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ABSTRACT 


Electron transport properties of terahertz (THz) longitudinal optical (LO)-phonon 
quantum cascade (QC) structures were modeled, in order to investigate high gain 
quantum cascade laser (QCL) structures. A new structure, a step well QC structure, was 
proposed. Under such an arrangement, there are three main energy levels within the step 
well, where the transition from the upper state to the middle state is at the THz radiative 
spacing and the transition from the middle state to the lower state is at or near the LO- 
phonon energy (~ 36 meV in GaAs). Because of the inherent difficulties in using rate 
equations for this type of transport analysis, a Monte Carlo simulation was developed. 
Step well injectors were modeled and shown to be capable of high injection efficiencies 
(~ 90%), higher than previously obtained. Comparisons to conventional square well LO- 
phonon structures are made, including a Monte Carlo analysis of a high power THz QCL. 
Interface roughness scattering was shown to be significant only for roughness greater 
than approximately one monolayer. It was found that step well structures are capable of 
high gains and injection efficiencies, with comparable characteristics to other square well 
designs, but do have increased scattering from the upper state to the lower states. 
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EXECUTIVE SUMMARY 


Electron transport properties of terahertz (THz) longitudinal optical (LO)-phonon 
quantum cascade (QC) structures were modeled, in order to investigate and better 
understand the underlying physics of high gain quantum cascade laser (QCL) structures. 
A new structure, a step well QC structure, was proposed after rate equation analysis 
showed favorable scattering rate lifetimes for maintaining a population inversion. Under 
such an arrangement, there are three main energy levels within the step well, where the 
transition from the upper state to the middle state is at the THz radiative spacing and the 
transition from the middle state to the lower state is at or near the LO-phonon energy 
(~ 36 meV in GaAs) for fast depopulation. The middle state (upper phonon or lower 
lasing state) is a single energy state, contrasting to previous LO-phonon based QCL 
designs that have doublet states. For long wavelength lasers, this may be important 
because it eliminates the possibility for unwanted THz absorption that could otherwise 
occur between those doublet states. Since the radiative and LO-phonon transitions are 
intrawell in nature, high oscillator strengths and sub-picosecond middle state lifetimes are 
possible, which can lead to increased gain in the active region provided the upper state 
lifetime and injection efficiency are maintained. Because of the inherent difficulties in 
using rate equations for this type of transport analysis, a Monte Carlo simulation was 
developed. Step well injectors were modeled and shown to be capable of high injection 
efficiencies (~ 90%), which is higher than previously obtained in other LO-phonon 
structures with similar wavelengths of operation and oscillator strengths. It is likely that 
step injectors could be useful in approaches that do not arrange all three energy levels 
within the same well. Comparisons to conventional square well LO-phonon structures are 
made, including a Monte Carlo analysis of a high power THz QCL. Interface roughness 
scattering was shown to be significant only for roughness greater than approximately one 
monolayer. It was found that step well structures are capable of high gains and oscillator 
strengths, high injection efficiencies, with comparable characteristics to other square well 
designs, but do have increased scattering from the upper state to the lower states. 
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I. INTRODUCTION 


Terahertz (THz) sources capable of operating in roughly the 1 to 5 THz region 
(300 to 60 //in) or far-infrared region are of considerable interest and are currently an 
active field of study. The THz frequency range has remained one of the least developed 
spectral regions, in part due to the lack of compact coherent sources. Though, in the last 
seven years or so, there has been a significant increase in efforts to produce compact THz 
sources. One approach to providing a compact source within the THz region, is to use a 
quantum cascade laser (QCL). 1 QCLs are solid state unipolar devices that differ from 
conventional semiconductor lasers in that transitions occur all within the conduction 
band, rather than between the bandgap of the conduction and valence bands. This class of 
semiconductor lasers allows for customization of the wavelength, by utilizing a multiple 
quantum well active region. A key feature of a QCL is that it consists of N repeated 
sections that form a cascade. Thus, one injected electron can emit many photons, which 
allows for differential quantum efficiencies greater than unity and hence higher power 
output. 

QCLs have been successfully designed for operating in the infrared (IR) region 
(for wavelengths from ~ 3 to 24 //m, with room temperature operation, and CW operation 
of hundreds of mW of power). - ' More recently, QCLs have been designed in the THz 
region (~ 0.84 to 5 THz, up to ~ 185 K, 250 mW pulsed, and 130 mW CW operation). 8 ' 12 
Though, for THz QCLs these figures are not from the same device. Temperature as a 
function of frequency for the best devices has approximately followed the trend T ~ 
Eradiativjk B = hcolk B , and several of the low frequency designs operate with the assistance 
of a magnetic field. Alternative compact solid state sources, such as lead-salt 
semiconductor diode lasers, are limited by the bandgap to ~ 30 jui n for the longest 
wavelengths. 14 Other solid state devices such as transistors, Gunn oscillators, and 
Schottky diodes multipliers, have at best achieved low //W power levels in the THz 
region. 15 ' 18 Gas lasers generally have limited lasing frequencies and are usually fairly 
bulky. For these reasons, and because of the lack of intersubband materials, QCLs 
provide a good means for solid state THz sources. 


1 



Some applications for the terahertz spectrum include astrophysical science, 
medical science, THz spectroscopy, various security applications such as THz imaging 
systems, chemical gas sensing, and agent detection. 15,19 ' 31 Because THz radiation is 
nonionizing, it is attractive for security imaging applications. Many materials have 
absorption bands in the THz region of the electromagnetic spectrum, therefore this may 
be used for the characterization of explosive materials. Imaging in this region can allow 
for characterization of materials as they are brought through a security screen, as well as 
detect and identify materials. THz waves penetrate dielectrics, such as clothes and plastic, 
but are blocked by metallic objects, which affect their screening and detection 
capabilities. 

A. PREVIOUS THz QC STRUCTURES 

THz QCL designs have used two main quantum cascade (QC) structure 

1 32 36 8 

architectures, broadly categorized as miniband ’ ' ’ and longitudinal optical (LO)- 
phonon, ’ ’ using coupled square quantum wells (square when unbiased). Although 
sometimes hybrid structures have been used with some overlap of these two approaches, 
usually designs can be classified as belonging to one of these two classes. 

In miniband designs, very closely spaced energy levels are formed which create 
minibands. Radiative transitions take place between these minibands or between a 
relatively isolated state and a miniband. The first THz QCL which used a chirped 
superlattice, 1 is an example of a miniband structure. Other miniband designs have 
included a bound to continuum (Figure I.A.l(a)) and a hybrid bound to continuum with 
optical-phonon scattering approaches. " Both featured minibands and somewhat isolated 
radiative states with more sideways radiative transitions. 

In THz LO-phonon designs, the approach is similar to the first QCL developed” 
(which was a mid-IR laser) in that the lower lasing state is depopulated by electron- 
phonon scattering. Though, no terahertz LO-phonon QCL was fabricated that used all 
diagonal transitions like the first mid-IR QCL. Also, the use of digitally graded alloy 
injectors has been replaced with easier to grow so called funnel or resonant tunneling 
injectors. ’ THz designs differ in that the radiative energy spacing is smaller than the 
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LO-phonon energy levels spacing. The first THz QCL that used the LO-phonon approach 
for depopulation, achieved the small radiative energy spacing by coupling the first two 
quantum wells with a thin barrier. This design approach still principally used three 
quantum wells, but the third well differed in that it was used to resonantly tunnel the 
lower lasing state (or equivalently the upper phonon state) and also arrange a lower 
ground state spaced near the LO-phonon energy (a fourth well was also included at the 
injector). Although some variations to these structures have been made, all of the THz 
LO-phonon based structures have essentially been within the same framework and have 
included this type of three well arrangement, differing only in the design of the injector. 




Figure I.A.l. Conduction band diagrams illustrating the two main broadly categorized 
QC structure architectures, with one section outlined, (a) miniband structure (showing a 
bound to continuum structure), and (b) LO-phonon structure (showing a three well 
structure). 
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This three well arrangement (Figure I.A. 1(b)) is exemplified in the given reference. 
Although at some THz frequencies, reported miniband designs have achieved the highest 
output power levels, this is attributable likely only to differences in the waveguide 
structure or the number of sections, and not the design approach of the active region. In 
fact, the highest output power QCL used a LO-phonon structure. 10 

LO-phonon based structures are of interest because of the large separation of the 
LO-phonon energy spacing, and also because the highest output power THz QCL used a 
LO-phonon approach. This approach seemingly is more likely tolerant to thermal 
backfilling, which could be important for increased operating temperatures. The LO- 
phonon structures typically have oscillator strengths (~ 0.5 to 0.96), while miniband 
structures such as the bound to continuum designs have higher listed oscillator strengths 
(~ 1.9). However, this may be somewhat offset since the LO-phonon sections are 
typically about half the length of their miniband counterparts, and thus typically have 
about twice the number of (LO-phonon based) sections within the same overall active 
region thickness. Our multiple section analysis of some miniband structures, has 
indicated lower oscillator strengths than reported. All of these miniband and LO-phonon 
designs have used square quantum wells that are symmetric when not under bias of an 
applied electric field. Next, a new theoretical approach that was investigated as part of 
this work, to the design of the active region for THz QCLs, will be discussed. 

B. STEP WELL QC STRUCTURES 

In addition to the use of miniband and square well LO-phonon approaches, the 
use of step well (asymmetric when unbiased) QC structures was first proposed by 
Freeman (author) and Karunasiri. 40 These types of structures can allow for the radiative 
and LO-phonon transitions to be placed within the same well. Since the bound state 
energy in high barrier square quantum wells increases ~ n" (with n being the quantum 
number of the state), it is not possible to have an upper radiative THz energy spacing and 
LO-phonon energy spacing (below) within a single square quantum well. Under a step 
well arrangement, there are three main energy levels, where the transition from the upper 
state to the middle state is at the THz radiative spacing and the transition from the middle 
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state to the lower state is at or near the LO-phonon energy (~ 36 meV in GaAs). The 
middle state (upper phonon or lower lasing state) is a single energy state, contrasting to 
previous LO-phonon based QCL designs that have doublet states. Allowing for a single 
middle energy state could be important for longer wavelength (lower THz frequency) 
devices to reduce unwanted absorption. By having vertical radiative and LO-phonon 
transitions within the same well, it is possible for these types of structures to yield high 
oscillator strengths, which can lead to increased gain in the active region provided the 
upper state lifetime and injection efficiency are maintained. The step in the well allows 
for high injection efficiency due to the spatial separation of the wavefunctions. 

A step quantum well, in which at least two different conduction band heights 
within a well are used, allows for additional freedom and breaks the restriction single 
square wells have. Thus, a THz and LO-phonon energy spacing can be arranged within a 
single step quantum well. Others have analyzed step wells for proposed C0 2 pumped 
THz laser applications. 41-42 Those step wells differ from the step wells considered in this 
research, where the radiative state is positioned above the LO-phonon transition states, 
and are intended for electrically pumped QCL structures. It was these unique 
characteristics of a step well, that prompted investigation of using a step well approach 
for a QC structure. 

C. OVERVIEW 

The ongoing imaging research using THz QCLs at NPS and by other groups, has 
shown there is a need for high power QCL sources. 43 " 50 The research documented in this 
dissertation focused on investigating high gain THz QC structures of a particular type, 
LO-phonon QC structures. This was because our step well structures, which are LO- 
phonon designs, showed characteristics favorable for a high gain active region, and 
because the previously referenced highest power THz QCL was a LO-phonon design as 
well. 

In the sections that follow, techniques developed for analyzing QC structures will 
be discussed. Methods for modeling quantum well structures, which are an absolute 
necessity for analyzing or designing QC band structures, will be covered. Various 
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scattering mechanisms and the calculation of pertinent scattering rates for detennining 
state lifetimes, will then be discussed. Fundamental rate equations will be also reviewed 
along the way. As will be shown, the scattering rates of the electronic states are not 
constant, but are density and temperature dependent. Because of this, rate equations are 
not well suited for analyzing complete structures, as they would be nonlinear. One 
approach that is suitable for modeling entire structures is one based on Monte Carlo 
simulations. The Monte Carlo method that was developed for modeling the electron 
transport in QC structures is then discussed. The last sections cover the design and 
analysis of step well QC structures, and the first sample that was grown and processed. 
Comparisons are made to conventional square well LO-phonon designs. The final section 
will show our Monte Carlo analysis of a square well high power THz LO-phonon QC 
structure, including the effects of interface roughness scattering, and comparisons are 
made to experimental data. 

The List of Symbols page shows some of the constants and notations used. Many 
of the derivations work in units with h = c = 1. However, all sections dealing with 
equations of final computable values will have h and c inserted back into the equations. 
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II. QUANTUM WELLS 


The energy states in the active region of a QCL are formed by quantum wells in 
the conduction band. In the sections that follow, methods for calculating the electronic 
states in quantum well structures are discussed. The transfer matrix method is developed 
for solving the energy levels and wavefimctions in such structures, both for arbitrary 
shaped potentials and for a linearly changing potential. The former is useful for analyzing 
doped structures, where bending of the band structure occurs due to the separation of the 
charges from the host ions. This requires a self consistent solution to Schrodinger’s and 
Poisson’s equations. These techniques are used to solve step quantum well structures, to 
illustrate their potential use for THz QC structures. 

A. TRANSFER MATRIX METHOD 

There are a number of ways to solve Schrodinger’s equation for bound states in 
quantum well structures. One commonly used approach, the transfer matrix method, can 
be used to solve arbitrary shaped conduction band potentials, which is useful when the 
self consistent solution is needed. Also to be discussed, is the case where the potential 
can be modeled as perfectly linear (classically a charge moving under constant 
acceleration). 

The most general solution can be formulated by considering the solution for the 
case where the potential is regionally constant. With this approach, any linear or 
nonlinear potential can be modeled by a succession of small steps. The solution to the 
Schrodinger equation H\y/) = E\ y/) for a region with constant potential can be written as 
a sum of forward and reverse waves (Figure II.A.l). 

y/{x) = Ae ikx + BAQ ikx (II.A.l) 

Since the wavefuncion must be continuous, at an interface (x n ) the first boundary 
condition can be written as 

¥ n i x n ) = y / n + lOJ (II.A.2) 
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Figure II. A. 1. Illustration representing the conduction band as a number of constant steps, 
showing the coefficients for use in the transfer matrix method. 


A second boundary condition, is that the first derivative must be continuous. Choosing to 
take into account the change in effective masses across the boundaries, it is noted that the 
conserved probability density flux is j = 1/m Im{ i//Vi//}. Setting/,, =j n +i the following is 
found 


1 dyf n {x n ) _ 1 ch// n+l (x n ) 

m n dx m n+ \ dx 


(II.A.3) 


The transfer matrix can be found from these boundary conditions by seeking relations for 
the coefficients A„ and B„ in terms of A n +1 and B n +\. Substituting equation (II.A. 1) into 
equations (II.A.2) and (II.A.3), the relationships for A„ and B n can be found. The resulting 
equations can be cast into matrix fonn 


A„ 


B„ 


1 + 


S,+1' 


k n m n+ 1 j 


1- 


K+\ m n 
k „ m n +1 J 


, i ( k , l + l - k ll ) x n 


AK+ i + i«K 


1 — 


n+V 


k„m n+ 1 ) 


1 + 


C +\ m n 


k„m 


.-'(i.ti+tiK 


-i(k„U~ k n) x n 


n"'n +1 J 


*n +1 


B, 


+i 


(II.A.4) 


2 1/2 

where inserting h back into the equations k„ = [2 m n (E-V n )/h \ . Because of the way the 
transfer matrix is developed, correct boundary conditions at the interfaces are always 
insured. In a structure with N regions, the matrix equation in terms of the first and last 
regions may be found as 


1 

o 

_1 

ii 

£ 

o 

£ 

i 

N> 

~ A N -ll 

- M 

Vil 


1 

3 

o 

o 

m oi" 

-1 

_1 

LV 





_ m io 

m„_ 

L^-.J 


(II.A.5) 


For bound state solutions, the wavefunctions must be square integrable and therefore 
vanish as x —> ±oo. From equation (II.A.5) we see the bound eigenstates can be found by 
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setting moo(£) = 0. The transmission coefficient can also be found for the case where an 
electron is tunneling through one or more barriers. Assuming an electron incident from 
the left means that Bn -i = 0. The transmission coefficient is then found from T(E) =j N -\/jo 
= k N -imo/(kom N -i)\l/moo\ . 

If the structure can be modeled as having a perfectly linear potential (as is the 
case for some QC structures where minimal bending of the band structure occurs from 
the doping), the transfer matrix can be formulated for the case of a linear potential. 
Writing the potential as V(x) = V n -e\E\x, Schrodinger’s equation becomes Airy’s equation 

d ¥n{X ] =P n (x) Vn (x) (II.A.6) 

dp MY 

2 2 1/3 

where p n (x) = (— |E| x-rj n )/j3, ij H = (E-V„)/e, ft= [h~ |E| /(2 m n e)\ , and the solutions are a 
linear combination of Airy functions. 

Wn W = A n Ai (Pn ( *)) + B n Bi (p n (*)) (II.A.7) 


The transfer matrix can be obtained in a similar manner as with the constant potential 
case, and can be simplified to 51 ' 52 


V B , 


111 

di(cc n+i )Bi \a n ) - —fy jM '(«„+ i)Bi(a„) 

w; + i 

m in 

Al (a„) A i '(«„ +1 ) - Ai '(a n )Ai(a n+l ) 

>K + i 


yyi 

Bi{a n+x )Bi —^73 Bi(a n )Bi '(«„ , ) 

™„ + i 

m in 

—jjy Ai(a n )Bi '(«„ +1 ) - Ai '( a n )Bi(a n+l ) 

w»‘+i 


B„ 


(II.A.8) 

where a„ = (2m„e\E\/h ) (-x„-rj„/ |E|) and the prime denotes a derivative. In a structure 
with N regions, we can again write the matrix equation in terms of the first and last region 
the same as in equation (II.A.5). For bound state solutions, the wavefunctions must again 
vanish as x —> ± 00 , or strictly speaking prior to the classical turning point as x tends 
toward 00 . Since Bi(a„) —* 00 as x —> - 00 , we set B 0 = 0, noting that mi 0 ~ 0 and mi 0 « 
mu. The bound eigenstates are found by setting mn(£) = 0. 

1. Self Consistent Solution of Schrodinger’s and Poisson’s Equations 
The self consistent solution of Schrodinger’s and Poisson’s equations is necessary 
to take into account doping in quantum well nanostructures. Formulating the transfer 
matrix method for a linear bias is computationally efficient, but is technically only valid 
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for undoped structures. In order to take into account bending of the band structure that 
occurs due to the separation of the dopant charges from the host ions, Schrodinger’s 
equation is solved by approximating the conduction band as a number of small steps 
using the transfer matrix method. This approach allows for solution of any arbitrarily 
shaped potential. Another popular alternative technique used is the shooting method, 
which deals with the effective mass at intermediate points by taking the mean at the 

53 

neighboring points. 

The solution to Poisson’s equation can be fonnulated for numerical computation 
as follows. Starting from Gauss’s law V-E = pi 8 and the electrostatic form of Faraday’s 
law VxE = 0, noting that the curl of the gradient always vanishes, the electric field can be 
written in the usual way as E = -V V a . This leads to Poisson’s equation V~V a = -pis. Now 
consider in general the following infinite sheet shown in Figure II.A.2, which is a plane 
surface of thickness Ax 3 with charge density mo(x), illustrating the electric field. 


side 2 


side 1 
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Figure II.A.2. Infinite plane surface of thickness Ax 3 and charge density « 3 d(x), showing 
the electric field. 


Applying Gauss’s law we find 
f d 2 .rE • n = — [ d 2 xp{x ) =— f d 2 xa(x) 

J r» J c J 


(II.A.9) 


where o is the 2D charge density. For this case, the electric field is then found from 
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(II.A. 10) 


a 
s 

The charge density as a function of the growth direction may be written as 
cr(x 3 ) = -e[n 2D t//* (x 3 )i//(x 3 ) - « 3 d (x 3 )]Ax 3 (II. A. 11) 

Using equations (II.A.10) and (II.A. 11), the electric field is found, and the potential can 
then be determined by integration 

V r7 =-jdx En (II. A. 12) 

If the dopants are concentrated in other states than just the ground state, the following 
modification is needed 

a(x 3 ) = (x 3 )<//(x 3 ) - n 3D ,.(x 3 )]lx 3 (II.A. 13) 

(=0 

where the sum of all the sheet densities for all of the states equals the total sheet density. 
The self consistent solution iteratively computes the perturbing potential and adds it to 
the conduction band profile until after a given number of iterations, a self consistent 
solution is reached. Equation (II.A. 13) assumes the charge distributions in the states are 
known. In practice, this is not known unless the electron transport of the full QC structure 
is analyzed and solved. One method for solving the state populations can be 
accomplished using Monte Carlo simulations, and this will be discussed later. 

For QC structures, doping is needed to introduce charge carriers within the active 
region, and also to help reduce the possibility of the formation of high field domains. 54 A 
typical LO-phonon QC structure may have on the order of ~ 200 cascaded sections, and 
an undoped active region would essentially represent a large intrinsic region. Minimal 
doping levels can result in lower threshold currents, and moderate doping levels have 
been shown to result in the highest operating temperatures. 55,56 In practice, only one well 
or a portion of a well or barrier is doped, in each section of the cascade. 

B. STEP QUANTUM WELLS 

As previously discussed, a step quantum well is not limited to have the bound 
state energy increase ~ rr as is the case for a conventional square quantum well. As such, 
the THz and LO-phonon energy spacing can be arranged within the same step well. This 
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approach is illustrated using the step quantum well structures shown in Figure II.B.l. 
These wells are comprised of Al x Gai_ x As layers with compositions of 
0.143/0.035/0/0.143. The step well thicknesses in mn are 29.9/12.3 for Figure II.B. 1(a) 
and 20.9/13.5 for Figure II.B.l(b). The energy difference between states 2 and 1 for 
Figures II.B.1(a) and (b) are E 2 \ = 13.7 meV (3.3 THz) and E 2 \ = 17.9 meV (4.3 THz) 
respectively. Both step wells have E l0 near the LO-phonon energy. The energy levels and 
wavefunctions were calculated assuming an applied field of 10.1 kV/cm (typical for LO- 
phonon QCLs) using the material parameters given in the references. 57 ' 65 This illustrates 
that it is possible to arrange different radiative transition energies and keep near resonant 
phonon transitions, within a single step well. This approach can be used to design 
radiative frequencies lower and higher than these as well. 



Figure II.B.l. Conduction band profiles of step quantum wells, (a) Conduction band 
profile of a step quantum well comprised of Al x Gai x As layers with compositions of 
0.143/0.035/0/0.143 and well thicknesses in nm of 29.9/12.3, with E 2 1 = 13.7 meV (3.3 
THz) and E i 0 = 37.9 meV. (b) Conduction band profile of a step quantum well comprised 
of Al x Gai_ x As layers with compositions of 0.143/0.035/0/0.143 and well thicknesses in 
nm of 20.9/13.5, with E 2 \ = 17.9 meV (4.3 THz) and £10 = 36.5 meV. For both (a) and 
(b), the 1 to 0 transitions are near the LO-phonon energy and the applied bias is 10.1 
kV/cm. 


This approach has a number of attractive features. Since the radiative and LO- 
phonon transitions are both vertical, this approach can yield large oscillator strengths and 
fast LO-phonon scattering rates for depopulation. Oscillator strengths for step well QC 
structures are typically around unity or greater. Conventional square well LO-phonon 
structures typically have oscillator strengths only about half of this and up to 0.96. It is 
also noted that the step represents an additional barrier that can reduce parasitic injection 
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into the middle state. This could help the injection efficiency in QC structures, especially 
for longer wavelength devices where the radiative energy transition spacing is small. It is 
also seen that the middle state is a single state, contrasting to other conventional square 
well THz LO-phonon designs that always have doublet middle states. Having a single 
state could help reduce unwanted THz absorption, which could be important for longer 
wavelength devices. In the following sections, the calculation of various scattering rates 
of the electronic states in quantum well structures will be discussed. 
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III. SCATTERING AND OPTICAL TRANSITION 


In order to determine the state lifetimes within QC structures, various scattering 
rates of those states must be computed. In the sections that follow, the calculation of 
scattering rates for the pertinent scattering mechanisms will be covered. As will be 
shown, the scattering rates and hence the corresponding lifetimes of the electronic states 
are not constant, but rather are density and temperature dependent. We will also discuss 
the optical transition, and the important calculation of the gain. Some approximate rate 
equations will also be used to illustrate fundamental device parameters of interest. 


A. ELECTRON-PHONON SCATTERING 


Phonons are quantized vibrations within a lattice. Many electrical and thermal 
properties of materials are governed by phonon interactions. In semiconductor 
heterostructures, electron-phonon scattering can often be the dominant scattering 
mechanism. If the spacing between two states is at or near the resonant LO-phonon 
energy, the scattering rate to the lower state can be very fast. This is why resonant 
phonon scattering has been used in LO-phonon QCL devices as the primary mechanism 
for keeping a population inversion. Intraband LO-phonon scattering along with electron- 
electron scattering, are the mechanisms that cool the electron gas within subbands and 
thermally distribute the electrons into Fermi-Dirac distributions. Since our focus is on 
LO-phonon QC structures, this is the most important scattering mechanism, and it will be 
discussed first. 

Beginning by considering a simple classical diatomic system, which is a good 
model for GaAs, the model consists of masses attached via springs in a chain (Figure 
III.A. I). 66 The Hamiltonian for the system is 

H = Y j ke + v 







(III.A. 1) 
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Figure III.A.l. Linear chain of diatomic atoms of masses nzi and m 2 . 


From Hamilton’s equations we can write 
8H 


P = — 


8x 


m x Xj = - 

nux, = - 

z J 


k(x' j -x J )+k(x j -x' j _ J] 
k(x'j-x)-k(x j+ i-x')] 


(III.A.2) 


or 

m 


l x j =k(xj+x' j _ l - 2 Xj ) 


m 2 x'j =k(x J+l +x j -2x' j ) 


J 

•r 

K j -\-- j +1 ■ "7 “" 7 / 

We seek solutions of the form 

1 ( 1 J a-wt) 


(III.A.3) 


x ; = RefA e 


x, = Ret A' e 


rt i {qia-at 


(III.A.4) 


where a is the lattice spacing between like atoms. Using equations (III.A.4), equations 
(III.A.3) can be written in matrix notation after some manipulation as follows (Re{} is 
assumed) 


(col -CO 2 ) 

< 

(l + Q~ iqa ) 

~X~ 

\ m 1 / 

_<(l + e -) 

2 

(or - 

\m, 

co 2 ) 

X' 


= 0 


(III.A.5) 


where 


col = 


2k 


m. 


col 


2 k 


(III.A.6) 
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and q = q—qj-x■ The solution is found by setting the determinant to zero. Solving the 

2 

quadratic in terms of co leads to the following. 


co 1 = 


(co 2 + CO 1 

\ m { m 


- + 
1 


f 2 , 2 A 2 

CO,„ + CO,,, 


2 2-2 

- co,„ m sin 


^ qa^ 


V z J 


(III.A.7) 


Equation (III.A.7) can be used to graph the dispersion curve of co versus q. Shown in 
Figure III.A.2 is the dispersion curve with m 2 = 0.93 m i (zwca = 69.72 g/mol, /«as = 74.92 
g/mol). As can be seen, there are two main branches of phonons, known as acoustic and 
optical phonons (the notation co —> coac for the acoustic branch and co —> coop for the 
optical branch, will be used). 



Figure III.A.2. Phonon dispersion graph, showing the acoustic and optical branches. 


For the acoustic branch, which corresponds to the “negative sign” root in equation 


1/2 

(III.A.7), as q —> 0, co —>■ 0 and we see using sin(ga/2) ~ qa/2 and (1 + stuff) ~ 1/2 (stuff) 


co... co „ 


®AC ~ 


qa 


7 


(III.A.8) 


+ a). 


Substituting equation (III.A.8) into equation (III.A.5) gives the relative displacement 
between the two atoms. 


Re 


\f + iqa^ 

z] 

ll 2 ) 

J 


(III.A.9) 
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For the optical branch, which corresponds to the “positive sign” root in equation 
(III.A.7), as q —> 0, a> —> nw, and when q = 0 the following is found 



(III.A. 10) 


where 1 lm re ducded = l/mi+l/»i 2 . Substituting equation (III.A. 10) into equation (III.A.5) 
gives the relative displacement between the two atoms. 


Rej-—f 1 + -^-1x1 = X' (IILA.il) 

[ m 2 V 2 J J 

As expected because of the high frequency nature at q = 0, the optical branch is 
characterized by opposing motion of the opposite atoms. 

From equation (III. A.7) it is seen that there are a range of solutions specified by q 
= ±7i/a. For periodic boundary conditions (x,+n = Xj and x]+n = xj) of a finite chain of 
length L = Na, i.e., the ends of the chain are tied together in a ring, qL = 2Tin, or q n = 
27vaJL. There are N number of m.\ atoms = N number of m 2 atoms = N number of two- 
body pairs or 2 N number of atoms total (which is the number of degrees of freedom for 
the system). For the 3D case, in addition to longitudinal phonons, there can also be 
transverse phonons. In mixed composition crystals, the main branches can split into sub¬ 
branches. This has been experimentally measured and theoretically derived using the 
Lyddane-Sachs-Teller (LST) splitting. 67 ' 69 

In polar semiconductors, the phonons create polarizations and thus electrostatic 
fields that follow the phonon. The interaction that describes this is referred to as the 
Frohlich interaction. An expression for this interaction can be found by first considering 
in general the dipole moment pi = q c x, where q c is the charge and |x| is the displacement. 
The dipole moment per unit volume V, referred to as the electric polarization P, is then 
written as 

P = lfK p i> (III.A. 12) 

The induced polarization due to the deformation by a phonon can be written as (noting e 
is the deformation charge, i.e., the charge that is transferred, not the electron charge) 
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(III.A. 13) 


P = e*(X-X') 


'\N 

V 


The continuity equation dj v = 0, where f ’= (p, J), can be integrated with respect to time 
to give 

dp 


J dtV • J = -^dt 


81 


(III. A. 14) 


Since P = \dtj, the charge density (charge per unit volume) is p = -VP. For an optical 
phonon with wavevector q, this becomes 


p(\, t ) = Re] - e*V • (X - x')— \ = Re] -1 eiq ■ (X - X')— e ,q ' x_i 




V 


J l V v yj 

The dot product q-(X-X') projects only in the longitudinal direction, i.e., the q direction, 
thus only longitudinal optical (LO)-phonons can induce a charge. The induced 
electrostatic potential (O^) can now be found from Poisson’s equation 
V“® q = -pis, where s= £», and noting from the divergence theorem 

\d 2 xV 2 $> q = -fd 2 x^~ (III.A.16) 


(III.A.15) 


thus 

koo V (iq) 


J e“ d 2 u | = Rej—^(X - X')— J e u d 2 


(III.A.17) 


where u = zqx -ia> q t, du = iq-dx, and d 2 u = \iq\’d 2 x. Perfonning the integral and taking 


the real part, 71 the following is found 

® =-l —— I(X-X')(< " ,v 1 e N 1 

“ 2 e n V q 


e' qx - e ,q x 


)= ————(X -X')(-e' qx +e^' qx ) (III. A. 18) 

2 ^oo V g 


and the perturbing potential is now found from 


V_ 


e -phonon 






- e' qx + e~ ,q 


■) 


(III.A.19) 


q 2 s „ V q 

An expression for the deformation charge e due to a LO-phonon can be determined from 


classical arguments, and is found to be' 


66 
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e = 


mVcofsJ 


n 

N 

Uoo 

£St J 


(III.A.20) 


To solve for the amplitude X-X\ the system must be quantized. Second 
quantization could be used to promote everything directly in our Hamiltonian to operators 
using the canonical commutation relations, but there is an easier way. This standard 
procedure is to write the Hamiltonian in terms of generalized coordinates. 66 Using 
canonical coordinates and momenta (Q/ and P/), a set of normal coordinates (q q ) is found 
such that the potential is in diagonal fonn. The canonical commutation relations can then 
be used to quantize the system in a straight forward fashion. The Hamiltonian expressed 
in canonical coordinates is 

V 2 




' 2m, y 2 


where for the restricted harmonic oscillator case 

d 2 V 


(III.A.21) 


A = 

ij SQ'SQ 


rQ'Qj 


(III.A.22) 


noting that A y is real and symmetric (Ay = A /7 ) and [Q/,P/] = iSy, [Q/,Q/] = [P/, P/] = 0. 
Following Feynman 72 by scaling our canonical coordinates and momenta (to make the 
notation simpler), we look to find the normal coordinates. The scaling chosen is Qi = 
(m ,) 12 Q/, Pj = P/(m,-) 1/2 and U,j = 2Ay /(«,■ mf) 112 . The Hamiltonian then becomes 

»ho = jzk+jz^ae, (m.A.23) 

^ i ^ ij 

noting that Uy is real and symmetric (Uy = Uj,). A transform between the nonnal 
coordinates and the canonical coordinates can be written as 

% = Z C „0 (III.A.24) 

i 

where the transformation matrix C q , is orthogonal, Zc of fSi - 5 oji" Y, C f‘-v - S ij , and the 

i q 

inverse transform is 

0 = Z C <A (III.A.25) 

q 

The canonical momenta is likewise written as (where \q q ,pp[ = iS q p) 
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(III.A.26) 


P,=Z C ^ 


The eigenvalues o\ can now be found from 

Xc„c s (/„=®^ 


and thus 




(III.A.27) 


(III.A.28) 


Using equations (III.A.26) and (III.A.28), the Hamiltonian can now be written as a sum 
of independent harmonic oscillators. 






rs x q J 1 

Z q Z 

f 2 . 2 2) 


(III.A.29) 


These decoupled oscillators represent noninteracting phonons. The solution is of course 
well know for the hannonic oscillator, and it is now straight forward to quantize. Treating 
each classical hannonic oscillator as a quantum harmonic oscillator, p q and q q are 
promoted to operators, and the canonical commutation relations are imposed. 
q ., p = id 

1 r y 1 (III.A.30) 

Hi^j\=\pnpj\ ={) 


The operators p and q can be written in terms of ladder operators, noting that the 
canonical commutation relation [q h pj] = iSy is the same as [a q ,aq] = 1. 


q -lk k+<) 


(III.A.31) 


The Hamiltonian for the sum of the independent oscillators can now be written as 


H HO =Z^q[ a l a q + |] = Z®q[ a q a !“|j (III.A.32) 

and the Hamiltonian can be explicitly written for the 3D case as follows, keeping in mind 
that for the 3D case both longitudinal and transverse modes can exist 
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+ — 


q 5=1 


(III.A.33) 


where for each independent oscillator |i = 0) is an eigenstate of Hho with a zero-point 
eigenvalue of l Ao) (l . When the zero-point energy is set to zero, a qjS |i = 0) = 0 and the 
vacuum state |0) has an energy eigenvalue of E = 0. Using the commutators below 


H HO’ a q,s J 


^q a q,.s 


| HHO’ a q,.s] 


^q a q ,s 


(III.A.34) 


it is seen |n) = (a qjS ) n |0) and the Hamiltonian for each oscillator is diagonalized, which has 
a complete set of eigenvalues of (n+*A )a>= (n +V 2 )E L o (for n = 0, 1, ...). 

The background energy of the phonon modes is 224 <u q . For systems with many 
degrees of freedom, this c-number becomes large and infinite in extent for fields. Since 
only energy differences can physically be measured, this zero-point value is not 
important for our purposes here. However, it is noted though for fields that have infinite 
degrees of freedom, this shift in the zero-point could potentially be a problem. 

The statistics of phonon particles can be determined by considering the following 
two particle state apa q |0). Noting that a p and a q commute, the state is identical if the two 
particles are interchanged. Analogous to a classical oscillator that can be excited to an 
arbitrary number of high levels, mode q can have an arbitrary number of particles. Thus, 
phonons are Bosons and the average number of phonons in a mode can be found from 
Bose-Einstein statistics ( N) = 1/Z Tr Ac /,7/ . 

The mode displacement is now found by using the correspondence principle 
where the classical energy in the mode is the same as the quantum energy in the mode, 
which is N q E w (where N q is the average number of phonons in the mode and E L o is the 
quantized eigenstate energy). 


\mco w 2 \X -Xf N = N q E w (III.A.35) 

\X - X'\ = j 2A ^° (III.A.36) 

1 1 \Nma>J 

Substituting the expressions for e from equation (III.A.20) and the modal 
amplitude X—X' from equation (III.A.36) into equation (III.A. 19), the perturbing potential 


for the LO-phonon interaction can be written as 
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where q = |q|, which is the magnitude of the 3D momentum vector (includes the parallel 
in-plane and X 3 components). Separating this into emission and absorption terms, with N q 
—> Nq+l for the emission case, to include spontaneous and stimulated phonon emission 



j. 1 /o 

This can also be written in terms of ladder operators (a q |M/) = (N q + 1) \N q +l), a q |A^) = 
(N q ) m \N q -l)) 



Sometimes this is expressed in tenns of the nondimensional constant ap = 
e 2 /(8 nfi co)(2m colfi) m {\l 1 / e st ). 

The following diagram (Figure III.A.3) shows the indexing and notation used for 
electron-phonon scattering, where i is the initial electron state,/is the final electron state, 
and n is the subband index. 
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f ri 

i, n 


\\ k|] ±qy {emission!absorption) 


/\/\/v 

phonon 


Figure III.A.3. Feynman diagram for electron-phonon scattering. 


In this section some vectors have three components (q) whereas others have two parallel 
in-plane components. Parallel subscripts will be used for parallel in-plane vectors (k|| and 
qn) and the prime notation will indicate the final state. Now with the Frohlich interaction 
potential, the scattering rate can be computed. The 2D wavefunctions are Bloch 
functions. 

\i) = \ w.a.) (III.A.40) 


Since the interest is in the matrix element for intersubband transitions, the following 

474 


approximation can be used 

(v f \ v \ Vi ) = { u f a r \ v \ 11 i a i )~{ u / \ v \ 11 i 


ce„\ a f a i) + [ U f U i 


a f \V\ a, 

cell \f 


a , 


V\a, 


(III.A.41) 


where for intersubband transitions (a)a) = 0 and {u{\u!) ce n = 1. Working in the (x| 
representation, where {x\, x 2 . Ay) = (x, v, z) with xi, x 2 parallel to the growth direction X 3 , 
the wavefunctions can be written as 


x z = 


VT 


e' k|r V„(^) : 


sfl 


e ik|r V„(x 3 ) 


(III.A.42) 


This will be used in subsequent sections and the periodic zz, part of the Bloch function 
will be dropped. 

1. LO-Phonon Emission 

Considering the emission case first, where V = V _ , . The matrix element can 

^ e -phonon 

be written as 
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(III.A.43) 

For mode q, there is only one value of qy that makes the d x\ integral nonvanishing, and is 
such when the exponent is zero (otherwise the other values lead to the addition of random 
phase terms which add to zero). This gives a Kronecker delta function, which is the 
conservation of momentum as expected for the emission case. 


lf\ V \i) = -L 
X ' A\ 2V 


1 1 


V ^00 £ st , 


^N q+ \ J dx 3 -e"* 3 * 3 rl,(x 3 V„(x 3 )A^( k n - k' - q (| ) 


eE, 


= -i 


LO 


2 V 


1 1 


\ £ o° £ S t J 




+ 




-*03*3 


q 

YAxi)V n (x 3 ) 


q 


^( k ii-k;-qii) 


(III.A.44) 

To find the electron-phonon scattering rate, the matrix element is substituted into Fermi’s 
golden rule, summing over the final states. 7:1 Inserting Pi and c back into the equations 
(noting E l0 = hco L o) 


2 n 


w_ , =— Y(/|f|i) \s\ 

e -phonon ^ I I /1 | 


E, + - 


trk'w ^ 
2m 


+ ha> w 


E..+- 


h 2 kA 


2m 


2n e fiat 


f \ 1 ' 


LO 


fi 2 V 




K+i)Z 


J <fx 3 


-'<73^3 




q 


^( k n — kf. — q.i) 


X - E n - E k) 

(III.A.45) 

Converting the final states summation to an integral in k-space 75 and dropping the 
momentum delta function 


25 




• iq ii 



Figure III.A.4. Contour illustrated for closing around the lower pole, for evaluating the 
contour integral in equation (III.A.49). 
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“ Q-m(.x 3 -x' 3 ) 


J dq 


3 2,2 

d\ +<h 


= 2m res = -2m 


-i(x 3 -*J)(—1?||) 

-- = j^ e -?i(*3-*3) for(x3 _ x ')>0 (III.A.50) 


-^ii +^ii q\\ 

and similarly for negative X 3 -X 3 , the contour is closed in the upper half plane 


00 -iq-,(x->—x\) 

£ -q» \x 3 -x ' 3 1 r* / t \ r\ 

3 —j-— =— e for(x 3 -x 3 )<0 


J dq 


q 11 +^3 d \\ 

and it’s seen that 


(III.A.51) 


°2 n 

■ = — e * 3 * 3 for all (x 3 - x 3 ) 


J dq 


3 2,2 


(III.A.52) 


^ 2 D (q||) = — J dx 3 dx' 3 </ n , (x 3 V„ (x 3 )<//. (xl )(//„ (xl )e q '" 3 
#11 


.76 


*7|| “*"^3 *7|| 

The / 2 D(qn) integral can now be written as 

(III.A.53) 

If we wish to include a screening term, equation (III.A.49) with screening becomes 7 

00 

7 2D (q pdsc) = J dx 3 dx' 3 y/ n ,(x 3 )y/ n (x 3 )yf* n .(x' 3 )ii/ n (x' 3 )J dq 3 K1 \ \ -^ 5 — (III.A.54) 

-00 

and the contour integral becomes 

00 

J d c h 
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= n\ 
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fe + oT 2 fa 2 + ql ) 2 fe + qi T~ 


(III.A.55) 


- V «if 1*3 —-*31 


Finally, upon substitution the hoiq^qsc) integral becomes 
ho (q N ,q K ) = n J dx 3 dx 3 y/\ (x 3 )y/ n (x 3 )yf* n . (x 3 )y/ n (x 3 ) 


I * 3 - *31 qj 


(q;+q 2 J ‘ 2 2 fe + ^) 2 faf+^r 


+? * c l- C 3 _JC 3l 


(III.A.56) 


Now consider the energy delta function in equation (III.A.47), which conserves 
energy. The inner function of the delta function can be simplified as follows 
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(III.A.57) 


E n ,+E k ,+hco LO -E n -E k = 0 


hoj,o-{E n -E n ,)+E k ,-E k =0 
M + ^k. q ~E k = 0 where ft ft/ = ftft7 /0 - (£„ 


fl0) 'e+^ n ( k \\ - 2k \\<l\\ cos 6 + 7||)- E k = 0 


where a similar notation to the given reference has been used. 77 Note that $ is the angle 
between k|| and q . The expression may be rearranged in tenns of the quadratic. 

, , ft 2 , . n 2 2 n 

h(D e - k n q n cos e + — = 0 

m 2 m 


q: -2k, cosOq,, + 


2 ma> e 
ft 


(III.A.58) 


= 0 


The roots to the quadratic are 


%± 


= k„ cos G±Aku cos 9 0- 


2 mco[ 


(III.A.59) 


and the following is noted 


77 


= 2 v k \\ cos 2 6 - 


2 m ft/ 
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(III.A.60) 


Using S{f(x)) = J]Six~x,)/\f (x t ) |, where Xi are the roots of fix), the delta function 

i 

becomes 




f Pi 2 fi~ 

- q, 2 - k, cos 6q„ + fim'„ 

y 2 m 2 m 
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ft 2 cos 2 0 - 
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2m 
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n \%e~% 

The final expression for the emission case is 


(III.A.61) 
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e -phonon 




IK + X )\ de T ---T^|^||/ 2 D (q||) [s(q n - q ]]+e )+s(q l{ - q lhe )] 


m+e %-e 


(III.A.62) 


Now consider the different emission cases that are possible. There are three 
different emission cases that can arise, and are referred to as cases la, lb, and 2. 

LO-phonon emission case la: Ek < fico' e = hco L o-(E ,- E„) and tico' e > 0, i.e., E n - 
E n ' < ha>w and the electron does not have sufficient in-plane kinetic energy, the transition 
is forbidden. The subbands are spaced such that the energy difference is smaller than the 
LO-phonon energy and the electron does not have enough in-plane kinetic energy, so the 
transition is forbidden (Figure III.A.5). 





Hcolo 



Figure III.A.5. LO-phonon emission case la, E n -E n ' < ha>w and the electron does not 
have sufficient in-plane kinetic energy, forbidden transition. 
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LO-phonon emission case lb: E k > hco' e = hco LO -(E n -E n ) and co' e > 0, i.e., E n -E n ' 
< ha>w and the electron has sufficient in-plane kinetic energy, the transition is allowed. 
The subbands are spaced such that the energy difference is smaller than the LO-phonon 
energy and the electron has enough in-plane kinetic energy, the transition is allowed 
assuming no state blocking (Figure III.A.6). The red and blue circles on the E-k diagram 
are equal energy planes and show the possible range for the in-plane momentum vectors. 
The final state represented by the blue circle must be vacant, otherwise state blocking 
will prevent such a transition. There are two possible solutions for this case, q\\± e . Limits 
on 0 are necessary to keep the momentum vectors real valued (to keep from getting a 
negative sign under the square root). Also note that at resonance the transition is allowed. 



Figure III.A.6. LO-phonon emission case lb, E n -E n ' < Hcolo and the electron has 
sufficient in-plane kinetic energy, allowed transition. 

Equation (III.A.62) becomes 
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(III.A.64) 
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LO-phonon emission case 2: Any Ek is valid and <x>' e < 0, i.e., E n -E n ' > ha>w, the 
transition is allowed. The subbands are spaced such that the energy difference is larger 
than the LO-phonon energy, and the transition is always allowed assuming no state 
blocking (Figure III.A.7). This “sideways” transition has one solution q\\+ e . When 6 is 
zero (Figure III.A.8) it is seen q = q\\+ e is the correct root, as q\\ must be greater than k . 


ha>Lo 




Figure III.A.7. LO-phonon emission case 2, E n -E n <> Pkdlo, “sideways” transition always 
allowed. 
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Figure III.A.8. LO-phonon emission case 2, qy and ky vectors when 6= 0. 


Equation (III.A.62) becomes 
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(III.A.65) 


2 . LO-Phonon Absorption 

The absorption case where V = F_ , , proceeds in a similar manner as for the 

1 e —phonon 1 

emission case. The main differences are N q +1 —> N q and the delta function changes for 
the conservation of energy for the absorption process. With these changes, similar to 
equation (III.A.47), the equation for the electron-phonon scattering rate becomes 
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e co 


LO 


e -phonon 


Stt 2 


1 1 
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N q \d6dq^ q ]} ^ 2D (q n ) S{E b , +E k , -hto w -E n -E k ) 

(III.A.66) 


Now the inner function of the delta function can be simplified as follows. 


E n . + E k , fi co lo E n E k — 0 


-tico w - (E n - E n ,) + E k , -E k = 0 

~k(o' a +E - E k = 0 where haf tt =ha> w - (E n , - E n ) 


(III.A.67) 


fi 2 / 7 X 

~h(o' a + —{k 2 - 2k n q l{ cos6 + q 2 )-E k = 0 


Again, 0 is the angle between k and q . The expression may be rearranged in terms of 
the quadratic. 


q , 2 -2k,, cosdq,, - 


h 


= 0 


(III.A.68) 


The roots to the quadratic are 
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(III.A.69) 



E„'-E n > Hcolo and the electron does not have sufficient in-plane kinetic energy, the 
transition is forbidden (for single phonon absorption). The subbands are spaced such that 
the energy difference is larger than the LO-phonon energy, and the electron does not have 
enough in-plane kinetic energy, so the transition is forbidden (Figure III.A.9). 
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Figure III.A.9. LO-phonon absorption case la, E n '-E n > ha>w and the electron does not 
have sufficient in-plane kinetic energy, forbidden transition. 

LO-phonon absorption case lb: E k > - hco' a = -[ ha>w-(E n ’-E n )] and aj a < 0, i.e., 
E n '-E n > ha>w and the electron has sufficient in-plane kinetic energy, the transition is 
allowed. The subbands are spaced such that the energy difference is larger than the LO- 
phonon energy and the electron has enough in-plane kinetic energy, so the transition is 
allowed assuming no state blocking (Figure III.A. 10). There are two possible solutions 
for this case, q\\± a . Limits on 6 are necessary to keep the momentum vectors real valued 
(from getting a negative sign under the square root). Also note that at resonance the 
transition is allowed. 



(III.A.73) 
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Figure III.A. 10. LO-phonon absorption case lb, E n —E n > ha>w and the electron has 
sufficient in-plane kinetic energy, allowed transition. 


Equation (III.A.72) becomes 
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(III.A.74) 
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LO-phonon absorption case 2: Any Ek is valid and <x>' a > 0, i.e., E n —E n < fra>w, 
the transition is allowed. The subbands are spaced such that the energy difference is 
smaller than the phonon energy so the transition is always allowed assuming no state 
blocking (Figure III.A. 11). This “sideways” transition has one solution q\\+ a which again 
can be seen when 0 is zero. 
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Figure III.A. 11. LO-phonon absorption case 2, E„<-E n < ficow, “sideways” transition 
always allowed. 


Equation (III.A.72) becomes 
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(III.A.75) 


LO-phonon absorption case 3: Any Ek is valid and (»[, > 0, and E n —E n < 0, the 
transition is allowed. Absorption from an upper band edge to a lower band edge is always 
allowed assuming no state blocking (Figure III.A. 12). This “upwards” transition has one 
solution q\\+ a . Phonon absorption scattering rates (all cases) are much slower than phonon 
emission scattering rates at low temperatures. This is because a phonon must be present 
to be absorbed. Further for case 3, absorption via a large momentum transfer also reduces 
the scattering rate. 
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Equation (III.A.72) again takes the same form as equation (III.A.75). 

3. Mean Scattering Rate 

Because there is a distribution of carrier energies in the initial state, the mean 
scattering rate is calculated by averaging over the Fermi-Dirac distribution of carriers in 
the initial state 


W , 

e -phonon, mean 



(III.A.76) 


where state blocking has also been included, and the electron-phonon scattering rate for 
both emission and absorption can be written succinctly as 
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(III.A.77) 


The populations and electron temperatures can vary greatly between the subbands, and 
these can only be determined by analyzing the entire QC structure. Modeling the entire 
structure using Monte Carlo simulations will be discussed later. 
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B. ELECTRON-ELECTRON SCATTERING 

Carrier-carrier scattering is a mechanism that can be significant, particularly for 
closely spaced energy states. When the energy spacing between states is smaller than the 
LO-phonon energy, LO-phonon scattering is forbidden except for electrons that have 
sufficiently high in-plane kinetic energy. Because of this, the mean electron-phonon 
scattering rate will be reduced. Electron-electron scattering, as well as other single 
electron scattering mechanisms (such as impurity and interface roughness scattering), 
may then be the dominant scattering mechanisms. Intraband electron-electron scattering 
along with LO-phonon scattering, must be modeled in order to determine the thermal 
electron distributions within subbands. For these reasons, modeling electron-electron 
scattering is important. 

The following diagram (Figure III.B.l) shows the indexing used for electron- 
electron scattering, where the electrons can be right or left handed, and i and j are the 
initial states and/and g are the final states. 


/ 


g 

j 


Figure III.B.l. Feynman diagram for electron-electron scattering. 



Spin independent interactions combined with the Pauli exclusion principle lead to a tenn 
in the form of SrS 2 , because S 2 = (S 1 +S 2) 2 = Si 2 +2SrS2+S2 2 . The spatial part of the 
wavefunction can be either symmetric or antisymmetric. 
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(III.B.l) 


The total wavefunction is the product of the spinor that describe the possible spin states 
of the two electrons and the spatial part. Since the electrons are fennions, the total 
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wavefunction must be antisymmetric under simultaneous exchange of both the spin and 
spatial parts. Using Pauli spin matrices 

nh 

(III.B.2) 


(III.B.3) 
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the spin can be written in the usual way as 
1 


S = —o 
2 


The permutation operator Pn = p n ip> ‘"p ]2 (spacc) has eigenvalue of 1 for the symmetric 
case and -1 for the antisymmetric case. For a two electron system, the eigenvalue for the 
SrS 2 triplet symmetric states is % and for the singlet antisymmetric state —A . This allows 
the permutation operator for the spin to be written as 72 ' 78 


U‘""=f(l + 4S,-S 2 ) 


(III.B.4) 


which can be rearranged to get 
2P n (spm) -1 = 4Sj • S 2 
Since 


(III.B.5) 


^1 '^2 ^ ®1 °2 


(III.B.6) 


the Pauli matrices are related to the exchange operator by (subscripts are the particle 
numbers) 

<V<t 2 =2P n (spin) -\ (III.B.7) 

The total spin and square of the total spin is 

S = ia 1+ i<T 2 (III.B.8) 

S-S = | + ^ 0 1 <U (III.B.9) 

From (III.B.9) and (III.B.7) the following is found 

S.S = l + ^ + ^o, • o 2 = P n (spm) + 1 (III.B.10) 
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The total spin for the symmetric case is 1 and for the antisymmetric case 0 (the 
eigenvalues of CTrCT 2 are 1 for the symmetric case and -3 for the anitsymmetric case). The 
total wavefunction may be written as 

K) = |£Vs) and Wa) = \^¥a) (IITB.il) 

Electrons are indistinguishable particles. Denoting the matrix element as M, the 
following subscript notation is used. 


M n =M^ 

=M j ^ f 


(III.B. 12) 


Antiparallel spin electrons do not interfere while parallel spin electrons do interfere. This 
is because Pauli’s exclusion principle states that no two electrons with the same quantum 
numbers (including spin), can occupy the same space. Since antiparallel spin electrons do 
not change sign upon exchange, the square of matrix element for the antiparallel spin 
case can be written as 

M 2 a =\M n f+\M 2l \ 2 (III.B. 13) 

Parallel spin electrons interfere and the wavefunction changes sign upon exchange. Thus, 

7Q 

the square of the matrix element for the parallel spin case may be written as 
M]=\M n -M 21 1 2 = M 2 -2M n M n (III.B. 14) 

where the minus sign comes about from the exchange. The matrix element is the same as 
that for the antiparallel case, except for the additional exchange term (negative signed 
quantity). This exchange effect results in a lowering of the scattering rate for parallel spin 
electrons as compared to the antiparallel case. 

If we assume equal distributions of antiparallel and parallel spins, the matrix 
element may first be written by including the 14 probability factors as 

M' 2 =±M 2 +^M 2 p =M 2 -M u M 21 (III.B. 15) 

Now since = Af,y^g/, an additional 14 factor must be included to keep from double 
counting. 

M 2 = = |( Ml-M, 2 M„)=^(M„M n +M 2I M 2 ,-M, 2 M 2 ,) (m.B.16) 
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1. Antiparallel Spins 

The first case considered is where the electrons are assumed to have antiparallel 
spins, so no interference occurs and the fonn of equation (III.B. 13) applies. All of the 
vectors (k vectors) in this section represent parallel component vectors (kj|) and the 
parallel subscript has been dropped to simplify the notation. The perturbing potential for 
the electron-electron scattering interaction is the Coulombic potential V = Vcouiomb- 
Working in the Born approximation and in the (x| representation, the following matrix 
element is computed where (jci, X2, X3) = ( x , y, z) with xi, X 2 parallel to the growth 
direction X 3 . 


fgV\i,j) = (fg 


a 

~sA 2 


a 


£X-X 




J 2 X|| J 2 xfyx 3 dx\\/ f (x 3 )// (x 3 )^/* (x 3 )i//y (x 3 ) - 


a '[(k, -X|| +kj -x[ )-(k r -ly +k g -x[)] 


x-x 


(III.B. 17) 


This can be rewritten as follows, showing explicitly the X 3 and in-plane parts of the 
integral separately. 


a 

~eA 2 


f r 7 1 */ x / \ / 7 \f j2 t2 7 

J dx.dxjf/jlx-, )(//, (x,)i//„(x 3 )i// / (x3)J d X]] d x„ 


>[(kf -X|| +k j -xf| J-fky. -Xy +k g -x[| ] 

c 



(III.B.18) 


In order to simplify this expression, the in-plane part which is essentially a 2D Fourier 
transform can be evaluated first. This can be done as follows. The in-plane part can be 
written as 


^d^Xyd 


2 x' 6 


<[(kj-X|| +k j -xj)-(k,--xy +k g -xj)] 


X[| x|i ~ + (x 3 x 3 ) 2 


= J d 2 x^d 2 x |j 


/Ak-Xn +/Ak'-x[| 


X|j + (x 3 x 3 ) 


(III.B. 19) 


where Ak = k— k/- and Ak' = k k„. Using the following Bessel function identity 1 


80,81 



the expression becomes 


(III.B.20) 
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2 j2 i e 

” J X, 


J d 2 x n d 


>[(k f -X|| +k j -x[| )-(ky -X|| +k s -xj)] 


X.-Xn + (X3 - X3 ^ 


= {j 2 x ll j 2 x,; e , ' Ak - x " +, ' Ak '- XIi {jg I, 




Jr\(Qw k 


= J J d 2 x^d 2 x' l dq n e 


/Ak-(X|| -X|| )+/Ak-X|| +/Ak'-Xy —q^ |x 3 — x ' 3 


J o(^||| x ||- x [||) 


Now define p = X||-X|j, and note that ax\\ = \dx\ = |J(x||-xy)| = \dp\ since 
with respect to the |dx||| integration, thus J 2 X|j = \dp\ = d(f>pdp 

00 

= \\d 2 x^dq ^' (Ak+Ak )x e x ^^d(j)pdp J^(q^p) 

0 

00 

= J J J 2 x | jJg || e !(Ak+Ak ) x ' 1 e” ?l1 |jt3 ~ x ' 31 1 d<fpdpQ iAkpcosi * , J 0 (q n p) 

0 

82 

Using the following equations “ 

271 

J d(fi bkp c°s^' Hlp = In i m J m (Akp ) 

0 

with m = 0 this becomes 

2tt 

J d0G iAkpcos * = 2nJ Q (Akp) 

0 

and 


CO 1 

J pdpJ 0 (Akp)J 0 (q\\P) = —d(q l ~ A k) 


di 


the equation now becomes 


2 j2 / e 

X, 


j[(k, -xy +ky -xj )-(ky -X|| +k s -xf|)] 


[ d 2 x,,d 

00 00 

= J J J 2 X||J^||e /(Ak+Ak } x e qi1 ^ 3 x ^2 /rj pdp e iAk O 0 (M/?)J 0 (</||P) 

0 0 

00 X 

= J J d 2 x[ dqf^^'Q^-^lx—Siq^ - A k) 


d\ 


_ ^ c ~giih 

d 


13 1 <i 2 x'e 


/(Ak+Ak')-x[i 


~ x [i|) 


(III.B.21) 
l] is constant 

(III.B.22) 


(III.B.23) 

(III.B.24) 

(III.B.25) 


(III.B.26) 
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with q || = A k = |k—k/| from evaluating the integral with the delta function. Now from the 
d~x || integral which has finite limits of integration, Ak = Ak' is the only value such that the 
integral will be nonvanishing, which yields a delta function representing the conservation 
of momentum. 


2 ” j2 x' 6 


J d 2 x^d 


i[(k f -X|| +k y -xj )-(k 7 -X|| +k g -ij)] 


In 


|x N — x'J + (xj — x 3 ) 2 

, 2nA 


(III.B.27) 


-9II |*3 -*3 


A(k / +k j; -k i .-k / .) 


= —e' ?l " 3 3 U^(Ak + Ak') = 

% % 

This agrees with the given reference. 76 This expression can now be substituted back into 

the original matrix element to arrive at 

2 ncA 


■?ll *3"*3 


(f g\ v \ l ’j) = dx 3 dx y f ( x 3 V, ( x i)¥* g ( x ',)¥j( x i)? 

+k g -k. -k ; .) 

which can be rewritten in the usual notation as (noting a = e /4n) 

(f g\r\<- J > = 4 *(?„) 5 (k, + K ~ k, - k ; ) 

where 

A i jfg (q ,|) = J dx 3 dx' 3 i/ f (x 3 }// i (x 3 )i//* g (x' 3 ) y/j(x' 3 ) e ^" 1 * 3- * 31 


(III.B.28) 


(III.B.29) 


(III.B.30) 


and again q\ = Ak = | k,—k/|. 

The scattering rate can be found by substituting the above matrix element into 
Fenni’s golden rule summing over the final states/and g. Inserting Pi and c back into the 
equations 


k-, =f a 


fg 


h 2 kl h 2 k 2 
E f + - L + E„+- s 


f Pl-k ^ 


( Pi 2 k p 


E.+ 

— 

E.+ 1 


{ 2m j 


l ^ 1,1 J 



2 n 
Pi 


fe 2 ) 

2 

y 

II ) 

U sAj 

fg 

% 


2m 2m 

A(k/ + k„ - k ; -k j)S(E Tf + E Tg -E Ti - E Tj ) 


(III.B.31) 

where E T is total energy, i.e., the subband edge plus the in-plane kinetic energy. 
Converting the final states summations to integrals in k-space 
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w 


ne 


A A 


his 2 A 2 (2 k) 2 (2ft) 2 
x 8(E Tf + E Tg - E t . - E Tj ) 


\d 2 k f d 2 k g 


Ajfgidw) 


% 


^( k / +k s _k <-- k /) 


(III.B.32) 




-> T I W 1 k,d 2 k„ 

32 ft 3 hs 2 J f 8 


Ajfgidu) 


<h 


S(k f + k g -k, -k j)8(E tj + E Tg -E Ti -E Tj ) 


Integrating over all the states of the j index electron and including Fermi-Dirac functions 
to take into account the occupancy distributions 


W- - =— 
e “ e 32 ft 2 h£ 2 


J d 2 k f d 2 k g d 2 kj 


A* to ii) 


% 


I 1 - f FD J ( k f)\ I 1 - fFD, g ( k g)\f F DJ % ) 


X S(k f +k g - k, -k y )8(E rf +E Tg -E T . -E Tj ) 

(III.B.33) 

Grouping the occupancy distribution functions into one term Fjf g , and evaluating the d 2 k g 
integral to get rid of the momentum delta function 

?n)f 


w 




x 8 


32 ft 2 fis 

r 4 - 2 7-2 h 2 k 2 

8 2 m 


d\ 


F 


jfg 


trk 

E f + 

2m 


f +E+'-^--E l - 


trkj E kE? 

2 m J 2m 


^- %:—j [ d 2 k f d 2 k j 

n 2 32 ft 2 hs 2] f 1 


AM 


% 




me 


7 \d 2 k f d 2 k ] 


AM 


F * s 


k} + kl-kf-q+^l(E, + E t -E,-E) 


k) + k]-kl-k) +| ”(e, + e ! -e i -e j ) 


16 ft 3 h 3 e 2j J J \ q {l 

(III.B.34) 

This equation gives the scattering rate of an electron with momentum vector k, averaged 
over the other initial particle states k ; . Since k g can be written in terms of the other three 
momentum vectors, the only unknown is k/ and this can be determined as follows by 
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using relative momentum vectors ’ 

k, = k,-k, (III.B.35) 

k fg = k s k f (III.B.36) 
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In this scheme, we look to replace the integration over d 2 kf with an integral over d 2 kf g . 
The following can now be written for k/ 


k /= k ,+ k ;- k ,= k ,+ k ;- k /,- k / 


-(k,.+k.-k /g ) 


(III.B.37) 


Since k f g is the only one of the two relative vectors that is a function of k f 


dk f .. =-—dk 


V.*, 


fg.X 1 


1 


dk f,x 2 =~2 dk f8,x> 

and thus 83,84 
d 2 k f = dk 


V 


/ . X\ dk f 5 x 2 


dk 


fg.x 1 


-dk 


A 


fg.x 2 


= -d 2 k 


fg 


(III.B.38) 
(III.B.39) 

(III.B.40) 


83 2 

as originally shown in the given reference. Now d~kf g = d6kfxlkf g , where 6 is defined in 
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Figure III.B.2, using notation similar to the given reference. 

It is often reasonable to approximate that the occupancy of the final states/and g 
are small, and ignore state blocking. This effectively means that F t f g ~ f FD j(kj [) since f FD / 
= f F D, g ~ 0. Including state blocking decreases the scattering rate only if the densities in 
the final states are appreciable. In order to solve for q\\, the law of cosines is used as 
follows 


u = ^+ k ) - 2k i k j cos (^ - a \) 
= k 2 + k 2 + 2kjkj cos(a[) 

k )um y = k f + k l ~ 2k f k g cos(x - P) 
= k) + k 2 + 2 k f k g cos(/?) 


and 

k 2 = k 2 + k 2 - Ikjkj cosla,) 

k l = k f + k l ~ 2k f k g cos (P) 


(III.B.41) 

(III.B.42) 

(III.B.43) 

(III.B.44) 
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Figure III.B.2. Momentum vectors used in the calculation of the electron-electron 
scattering rate. 


Summing equations (III.B.41) and (III.B.43) and summing equations (III.B.42) and 
(III. B. 44) 

2 , ,.2 _ 0 7,2 , o ,,2 (HI.B.45) 


Kumij + Kj = 2kf + 2kj 


+2 K 


Substituting equations (III.B.44) and (III.B.45) to get rid of k 2 mm 
k 2 fg = kfj+2(k 2 f +k 2 g -kf-k>) 


(III.B.46) 


(III.B.47) 


The energy delta function yields conservation of energy. 

-2 K 2 

4- " ’ Vg = F+- 
"g 


E f + 


fi 2 kl h 2 k 2 n 2 k 2 n 2 k 2 

j _e -\ _£. = E + 1 + E +■ 1 


2m 


k}+k;-k 2 -k 2 = 


2m 

2m 


2m 2m 

(e i+ e,-e,-e s ) 


Substituting equation (III.B.49) into (III.B.47) 


(III.B.48) 
(III.B.49) 
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(III.B.50) 


k 


2 

fg 


= kl + 


Am 


(e 1 +e i -e ; 



+ k 


2 

O 


2 2 

using the usual definition k 0 = 4m/h (E,+E—Ef-E g ). Recall from the definitions of q\\, k tJ , 
and kf g , we may write 

|j^ _I 

=|k.-k / .|= 1 0 /g| (III.B.51) 

Now from the law of cosines 

(2<7h) 2 = ky + k 2 fg - Iky k fg cos 9 (III.B.52) 

Substituting kf g from equation (III.B.50) into equation (III.B.52) gives 

(2q/ = 2 k 2 +k 2 -Iky -\jky + k 2 cos0 (III.B.53) 


2 

Performing the dkf g part of the d kf g = ddkf g dkf g integration and using equation (III.B.47) 
gets rid of the delta function and gives q\\, which is a function of k/ g and 6. Thus, the d6 
part of this integration remains. 


W 


me 


me 


4\67T 3 h 3 £ 2 
\ 

jddd 2 k 


J d 2 k fg d 2 k j 




ufg 


4 \\ 




2 , ,2 ,2 , 2 2 m 


k~ f + k g - kr - kj+\E f + E g - E. 


( E f 


64 x 2 h 2 s 2 




% 


F 


jfg 



(III.B.54) 

Now d 2 kj = da\kjdkj, where a\ is the angle between k, and k ; (the integral is performed 
along kij with angle a\). 


W 


me 


e ~ e Mntis 2 


J dOda x k. dkj 


Aifgiy II ) 


% 


F 


jfg 


(III.B.55) 


The upper limit of integration for kj can be taken to be the top of the well barrier. 


k = \ 2m (v - F ) 

^ j max 'y ^2 x barrier ^ j ) 


(III.B.56) 


Equation (III.B.55) is an expression for calculating the electron-electron scattering rate. 

2. Screening 

The formulation so far has assumed that there are only the two initial electrons 
and they scatter due to the Coulomb force between them. That is, it has been assumed 
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that there are no other electrons (free carriers) present. In reality, there will be carrier 
concentrations in the subbands. These free carriers respond under the influence of the 
electrostatic field. This will have the effect of reducing the force between the initial 
electrons and reduce the scattering rate due to this screening. 

One screening model, is the static single subband screening model, where the 
dielectric function is modified to take the screening into account. This was first derived 
by the given reference and later used by others. This model can be included by taking 
s inside the scattering integral, and replacing it with s = ss sc , where 

+ (III.B.57) 

2s q 


and the polarization factor ignoring collision broadening at T— 0 K is 

\-\]{q,-2K F \\l- 


n„(a,r = 0 K) = 


m 


7itr 


1 


f \ 2 

2 K/ 


V 4\\ J 


(III.B.58) 


where \J(q\\-2KF) is a unit step function that when q\\> 2 K f is equal to 1, and otherwise 0. 
At any temperature T, under the assumption of noninteracting polarizability, the 
polarization factor becomes 




n„(«„r = 0K) 


= n„(<?„r = 0K)- 


4 k B T cosh' 


f E f -E A 

v 2 k B T J 


(III.B.59) 



f e f \ 

1 + tanh 

F 


V 2 k B T J_ 


where kn is Boltzman’s constant and Ki, is the Fermi wave vector defined as 


K f = 


2 n n. 


g 


(III.B.60) 


with the valley degeneracy factor g = 1 for GaAs, which gives the number of equivalent 
energy bands. 

As an example of how screening affects electron-electron scattering, consider the 
1100 scattering process for an infinite GaAs well 40 nm in length (Figure III.B.3). The 
rates are calculated for subband populations of 1x10 and 1x10 " cm at T = 300 K, 
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unscreened and including screening. As expected, the scattering rate increases as the 
subband population increases. The screening effect becomes significant for large sheet 

i o _o 

densities. Sheet densities on the order of 10 cm “ are two orders of magnitude larger 
than typically used in QCLs. 




Figure III.B.3. Electron-electron scattering rate as a function of total energy at T = 300 K, 
with and without including screening, for subband populations of (a) lx 10 10 and (b) 
lxlO 12 cnT 2 . 


3. Mean Scattering Rate 

The electron-electron scattering rate calculated so far, is for the i indexed electron 
at a given energy averaged over the j indexed electron initial state distribution. Since 
there is a distribution of carriers in the initial state, the mean scattering rate averaged over 
the Fermi-Dirac distribution of carriers in the initial state can be found as follows 


W 


e —e , mean 


since 


^ y ro ,(*, > m 




(III.B.61) 


>h = J dE k f FDi (k t ) 


m 

7th 2 


(III.B.62) 
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4. Antiparallel and Parallel Spins 

In this section, taking into account spin and state blocking will be discussed. 
Starting with equation (III.B.33), the solution can be found by defining angles between 

on 

the momentum vectors different than what was used earlier. Including screening, the 
equation can be written as 


W = —-- 

^ 32 7T 3 hs 2 


J d 2 k f d 2 k g d 2 kj 


A afM ii) 




r» 


[i -/™,/M 1 - V A >1 w*;) 


(III.B.63) 


x J(k/ +k g - k, -k j)S(E Tf + E,, g -E Ti - E rj ) 
rating over the d 2 k g to gc 
da\kjdkj, the equation becomes 


Integrating over the d 2 k g to get rid of the momentum delta function and using d 2 kj = 


W _ = 


- J d 2 k f da x k j dk j 




% £ s> 


[l -f FD A k M-f F »J^V FD j(kj) 


'bljr'ks 1 
xS(E Tf +E Tg -E T .-E Tj ) 

From the usual conservation of energy and momentum we have the following 


n 2 k) n 2 k 2 trk 1 ti 2 k 2 . 

E f + — f - + E+ -^ = E + —i- + E + — J - 

2m " 2m 2m 2m 


‘f 


k g =k i + kj-k f 


Now defining the angles for the momentum vectors as follows (Figure III.B.4) 


(III.B.64) 


(III.B.65) 
(III.B.66) 



Figure III.B.4. Momentum vectors used in the calculation of the electron-electron 
scattering rate, for the case when state blocking is included. 

2 

k g can be solved for by squaring the conservation of momentum equation 
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(III.B.67) 


K = ( k / + k 7 - k / )• ( k / + k , “ k / ) 

= k 2 + k 2 + k 2 f - 2 k t k f cos(y) + 2 kftj 005 ( 0 ;!)-2 kjk f cos(y -or,) 
Substituting equation (III.B.67) into equation (III.B.65) 

n 2 k 2 n 2 k 2 n 2 k 2 f n 2 

E. + —^- + E,. + — 3 - = E, + — j - + E+ — 

2m 2m 2m ‘ 2 m 

x \k 2 + k 2 + k 2 - 2kjk f cos(y) + 2k j k j cos^) - 2k j k f cos(y - or,)] 
The resulting quadratic can be solved to find kj 


(III.B.68) 


1 


1 


1 


\\k t cos(y) + kj coslcf! - y)f - 4k j k j coslc^)] 


k f = —kj cos (y) + - kj cos («! -y)±- 


4m 


(■ E,+Ej-E,-E ') 


(III.B.69) 


From Figure III.B.4 the following equations for k„ and q 2 are found, where angle y is the 
angle between k, and k/ 

k g = [(&,. + kj cos (a^-kj cos (y)f + (kj si n(a x )-k f sin(y)) 2 ]^ 


(III.B.70) 

= kj +kj-2kjk f cos(y) (III.B.71) 

These equations are in a form so that state blocking can be included. Now to take into 


account spin, consider the case where the electrons have equal probability of having 
parallel and antiparallel spins. Including the exchange effect discussed earlier, it is seen 
that the following substitutions must be made to equation (III.B.64) 

(III.B.72) 

?ii)| W*cW\\)\ WsM\\) £ sM\\)\ 

where 

q' 2 =k 2 +k 2 f - 2kjk f cos(«! - y) (III.B.73) 

and angle y-a\ is between k/and k ; . This agrees with that in the given reference. 88 


AM) 

2 

1 

_ ^ _ 

4,/g(<7||) 

2 

AM) 

AM^Am^) 

V\\ £ s M\) 

2 

% £ sM\^ 

I 


q\\if s M\) £ sM[) 


C. IMPURITY SCATTERING AND INTERFACE ROUGHNESS 
SCATTERING 

In this section, two additional scattering mechanisms, impurity scattering and 
interface roughness scattering, will be discussed. Beginning with impurity scattering, 
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ionized impurities within a semiconductor lattice and its distributed charges, can cause 
scattering due to the Coulombic interaction. The corresponding matrix element may be 
written as 


i/r !'■>=</ 


a 


£ X - X 


a 

~sA 2 


J d 2 x^d 1 x'^dx i \f/* f {x i ) \f/ i (x 3 )- 


a *(kj -k f )-x 


x-x 


(III.C.l) 


This matrix element can be simplified for scattering from an ionized impurity at X3 = X3, 
in a similar manner as was done with the electron-electron scattering case, inserting h and 
c back into the equation this becomes 

(/TO = (* 3 ) e ”^* 3- * 3 

CJl (III.C.2) 

e 2 4/(g||) 

2sAq n 

where 


4M\) = J dx 3 ^/(x 3 V/ (* 3 ) e~ ?lll * 3 “* 31 (III.C.3) 

and q\\ = A k = |k— k/| = k 2 + kj -2k t k f cos#. The scattering rate can be found by 
substituting the above matrix element into Fermi’s golden rule, summing over the final 
states 


K„=YlL\[rn)'A 


2 n 
h 


f e 2 * 


2sA 


f 


Z 


2 7,2 


E, + - 


trk 

_ j_ 

2m 


E : + 


h 2 k 2 ^ 

2 m 






(III.C.4) 


S(E Tff - E Tj) 


where E T is total energy, i.e., the subband edge plus the in-plane kinetic energy. 
Converting the final states summation to an integral in k-space, this becomes 
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\y = ^ ne _ d \d 2 k 8(E -E ) 

'« Me 2 A 2 (2x) 2 > f q, Cr/ T ’ f) 


8 7th s A 


\d 2 k f Af ^ S(E Tf - E Ti ) 

J J n J 


(III.C.5) 


4 jtfe z A 


J d 2 k f 


4 (?i|) / / 2 ,2 , 2 m 


8 k 2 -k 2 +- r {E f -E,) 


Performing the dkf part of the d 2 kf = dOkjdkf integration gets rid of the delta function and 
the following is found 


A , f (di) 


impurity o 3 2 


8tt/z ^ 4 • 


(III.C.6) 


To apply this for the case of doping over a distribution, where a sheet charge is dx^n(xi) 
at X 3 = JC 3 , the final expression for computing scattering due to ionized impurities becomes 


impurity 8 jtfe 2 


J dx[n m (x[)^dO 


A >M ii) 


(III.C.7) 


Monolayer fluctuations are often formed at heterostructure interfaces. Scattering 
can occur from these imperfect growth surfaces, where these variations in the barrier 
thicknesses gives way to variations in the energy levels and wavefunctions of the 
subbands. Assuming that the roughness height A(x) at the in-plane position has a 


correlation function of 


A(*,)A(*;))=Af( 1+ «A 


(III.C. 8 ) 


where the A is the roughness and A is the correlation length, the matrix element is given 
by 




(III.C.9) 


F mn = V 0 Vf(X 0 )Yi CO 


(III.C.IO) 


Substituting the matrix element into Fermi’s golden rule and summing over the final 


states we find 
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This is the scattering rate due to interface roughness using an exponential autocovariance 
function. A roughness a/2, where a is the lattice parameter, with a correlation length of 5 
nm was used in some of the simulations that will be discussed later. 

1. Mean Scattering Rate 

Because there is a distribution of carrier energies in the initial state, the mean 
scattering rate is again calculated by averaging over the Fermi-Dirac distribution of 
carriers in the initial state. 


l dE J v i„ 


impurity, 

interface roughness 


interface roughness 


ughnessj’ 

\dE ki f FDJ {k) 


(III.C. 13) 


D. RADIATIVE TRANSITION 

The optical transition occurs when an electron releases a photon and transitions 
from the upper to lower lasing state. As will be shown, spontaneous emission lifetimes 
are on the order of microseconds, while other rates previously discussed such as phonon 
rates can have lifetimes on the order of picoseconds. Thus, state lifetimes are not 
determined from spontaneous emission rates. This section will discuss the optical 
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transition and the very important parameter the gain. It will conclude by discussing some 
fundamental parameters found by approximate rate equation analysis. 

The potential for the radiative transition can be found by starting with the 

82 

nonrelativistic electromagnetic Lagrangian, and then finding the Hamiltonian in the 
usual way H = Y J P, t l, ~ L em- 

H = - 1 - (p 2 - 2qA • p + q] A 2 )+ q c O « + q c O - — A • p (III.D. 1) 

2m 2m m 

The Hamiltonian is now in the form of H = //<,+ V. Substituting for the electron charge, q c 

= -\e\, the interaction is given by 79 

V radiative =~ A P (III.D.2) 

m 

The equations for the E and B fields in terms of the scalar and vector potentials (<b and 
A) can be written as a second rank antisymmetric field strength tensor F MV = d^A v -d'A^. 82 
These equations automatically satisfy Maxwell’s homogeneous equations, and the 
inhomogeneous equations are d M F MV = J v . The scalar and vector potentials are not 
uniquely determined because the following change A M —*■ A fl +d fJ A(x,f) has no affect on 
the fields E and B since d M A v '-d'A M '= d M A v -d v A M . Therefore, one such gauge 
transformation we are free to choose is the Lorentz gauge d^A^ = 0, which Maxwell’s 
equations simplify to one equation UA M = J u . Without sources J = 0, and the equation for 
the vector potential becomes 
7)2 A 

V 2 A-^—= 0 (III.D.3) 

dt 1 

which has a solution of 


A = A o cos(q • x - cot)e = + e - i(q x - <w) ]e 

Since O is included in H 0 , E and H are solved effectively with 0 = 0. 

zv 

E =-= ~o)A o sin(q • x - cot)e 

8t 

1 A 

H = — Vx A =-s-sinCq-x-ruOqxe 

M F 

The Poynting vector can now be found from these fields 


(III.D.4) 


(III.D.5) 
(III.D. 6) 
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coA ~ 

S = ExH = —-sin 2 (q -x-cot)q 
M 

and the average Poynting vector is (noting co = v-q) 

S) = I^q = ^ v 


(III.D.7) 


. . (III.D.8) 

2 // 2 // 

Inserting fr and c back into the equations, the average energy density found from equation 
(II.D.8) is equal to Nhco/V, where V is the volume 


Ey) = 


S) A]q 2 Ntico 


v 2/u V 

from which A~ 0 can be solved for 

. 2 2 fiN 

A: = 


Y SCO 

Substituting this into equation (III.D.4), the vector potential becomes 
1 12 fiN\ 


A = —, 


k 


(q -TL-at) e -i(q x-frtf) g 


(III.D.9) 


(III. D. 10) 


(III.D.l 1) 


2 V V sco 

This can now be substituted into equation (III.D.2), the expression for the radiative 
potential. This is valid for absorption, and the following change N —* N +1 must be made 
for the emission case, to account for spontaneous and stimulated emission. In tenns of 
ladder operators, the interaction potential for the radiative transition becomes (dropping 
the harmonic part and promoting the momentum term to the quantum mechanical 
operator) 


F 


n 


radiative 


[a 1 e i '-+a;e- i '->-P = 


h 


m \ TV sea 


[ a q + a l q ]e iqx e-p (III.D.12) 


m V TV sco 

1. Emission Rate 

To find the 2D spontaneous emission rate, equation (III.D.12) is substituted into 
Fenni’s golden rule keeping the emission tenn. 




ne 


(III.D.l 3) 
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Using the electric dipole approximation, e iq x ~ 1 

+ l)Z|e ■ (/|p|')f <y ( £ > - E, - hm) (III.D.14) 

dVm co 1 

Setting N ~ 0 which is permissible for low intensity light, that is spontaneous emission 
occurs without any photons in the cavity (to stimulate emission) 

K: = ~ e < -M (III.D.15) 

sVm o)f' 1 

Summing over the polarizations and final states (including all final photon states for 
spontaneous emission into all modes), and converting the k-space summation into an 
integral 

w „..= ,, A ^ 'L\<l 2 p'fif\9\^S{E f -E l -hco) (III.D.16) 

dVmco {2k) h 11 

Note that d 2 p = dfipdp = d(f>'hoo/v 2 d(froo) 
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spo " 4sVm 2 KV 2 ti 


d</>d(ho))\e-S(E f -E. - tied) 
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4sL 3 m 7rv h pol 




(III.D.17) 


To simplify the bracket, for the moment we switch to the Heisenberg picture, where the 
operators (fix) = (/AxJ) = e"'(fix)e and dx/dt = \/(ifr)[xJ!] = p/m. 78 Substituting this 
commutator into the bracket 

(/| P |;) = ;^(/|[if,x]|;) (in.D.is) 

Noting that the momentum operator is a constant of motion, i.e., p(t) = /;(()), switching 
back to the Schrodinger picture we find 


(f\p\i) = ij{f\Hx-xH\i)=ij(E f -E,Xf\x\i) (III.D.19) 

and defining co = ( E—Ej)/ti (which will give a minus sign) 

(f\p\i) = -inuo(f\x\i) (III.D.20) 
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Substituting this into equation (III.D.17) and evaluating the summation for the 2D case, 
noting the polarization sum does not affect anything for the 2D case since the integration 
is only in terms of d(j)' (i.e., not d6) 


W =--- 

spo " 4sL,m 2 7iv 2 h 
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m cd 


pol 


2 2 

e co 


(III.D.21) 
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The final expression for the spontaneous emission into all modes becomes (noting 
envelope wavefunctions are assumed) 
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spon 
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e cd 
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(III.D.22) 


The 3D emission rate can be found in a similar method as for the 2D case, except when 
performing the polarization sum for the 3D case the square of the bracket tenn reduces to 


\(f\x\i)\ 2 sm 2 0 
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(III.D.23) 


As can be seen, there is a difference between the 2D and 3D emission cases (noting both 
are single electron scattering rates). A quick computation reveals spontaneous emission 
rates are on the order of microseconds, and thus as previously mentioned, do not affect 
state lifetimes. 


2. Gain 

The net transitions from the upper to lower lasing state, lead to an induced power 
= hco (N 2 Wi\-N\ W | 2 ) = hoj(AN x W 2 i), where W 21 = Wn are the stimulated emission and 
absorption rates respectively. This amplifies the propagating electromagnetic waves 
intensity as dl/dx = gain x l. The optical gain of the medium can be found, starting with 
equation (III.D.14) by keeping the stimulated emission term. Defining the photon flux 
Iphoton as the number of photons per unit area and time, substituting for the momentum 
operator and noting IphotoJv = (Number of photons)/V (where v is the speed of light in the 
medium and V is the volume), the stimulated emission rate (a positive quantity) assuming 
single mode becomes 
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(III.D.24) 


The gain is defined as (number of transitions per unit volume and time)/(emitted photon 
flux) = WILAtI photon = WlLAI p hot on . Noting that the population inversion is Aiv, n LA number 
of electrons, the expression for gain becomes 


. 7TC~fa Ipho 
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Taking into account a finite line width, the delta function can be replaced with the 


Lorentzian 


£(Hco) = 


FWHM/2 


x {AE-hoof + (FWHM /2) 2 


(III.D.26) 


where FWHM is the full width half maximum. The maximum gain occurs when AE = hoo, 
and the Lorentzian is then equal to 2!{n FWHM). Substituting this in place of the delta 
function, the expression for peak gain becomes 


, . 2 eh A«, n t r \ d 

peak gain = —;-—— ( / — 

m 2 sv6oFWHM V 1 dx 
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(iII.D.27) 


In Heaviside-Lorentz units sv = n inc i ex c (while in rationalized mksA units sv = s 0 n inc / ex c). 
This expression can be used for estimating the threshold population inversion, when the 
threshold gain is found from the waveguide resonator analysis. While having a large 
oscillator strength = 2h/mco\(J\d/dx\i)\ = 2mao/h\(J\x\i)\ is favorable, so is having a long 
upper state lifetime to maintain the 3D population inversion A» 3 d, as can be seen from the 


equation for peak gain. These are often competing characteristics. 

3. Rate Equations 

As we have seen, scattering rates are density and temperature dependent. Even to 
compute the mean rates, requires knowing the densities and electron temperatures of the 
states, which can only be detennined by analyzing the full structure. Thus, rate equations 
are not well suited for analyzing full structures. Nevertheless, some important 
fundamentals can be found from an approximate rate equation analysis. The rate 
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equations for a three level system, in terms of 3D volumetric densities, for the optical 
mode and state populations are 


(n 2 -n 1 )TVn p 


1 sport single mode 


T T 

spon all modes photon 


J (n 2 ~n,)Wn r 


1 spon single mode 


1 spon all modes 


dn i _^ ^ J + (« 2 ~«i)V% , 
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(III.D.28) 


(III.D.29) 


(III.D.30) 


1 spon single mode 


" spon all modes 


where the upper and lower lasting states are indexed by 2 and 1 respectively. The 
notation used is nphoton is the 3D photon population mode density, r spon sing i e mode is the 
spontaneous emission into the single lasing mode, T spon ail modes is the spontaneous 
emission into all modes (noting 1 / z spon a it modes = N modes 1 lr spon sulg i e mode and the total 
stimulated emission rate 1 l^stim 'Z n photott 1 /Lpon single mode)i r'photon r( (iro (hot) is the 
photon cavity lifetime, r m = va m is the mirror escape time, N is the number of sections, V 
= LA is the volume of one section, V c is the effective volume of the cavity, T = NV/V C is 
the confinement factor, 77 is the injection efficiency, and J = I/A is the 2D current 
density. At steady state with r 2 « r spon a ii modes (since r 2 i « T spon a ii modes) these 
equations become 


Q _ ( t h- f h) T ^ n photon { 


1 spon single mode 


T T 

spon all modes photon 


J (n 2 -n x )Vn hoton 
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spon single mode 


0 = ( 1 - ? 7W + 


J | (n 2 -n l )Vn photan | n 2 n x 


(III.D.31) 


(III.D.32) 


(III.D.33) 


1 spon single mode 


Above threshold, lots of photons are present in the cavity and Vn p i 10 ton is large, 
and the second term on the right hand side of equation (III.D.31) can be ignored. 
(n 2 -n,)rv _ 1 


T T 

spon single mode photon 


(III.D.34) 
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From equations (III.D.32) and (III.D.33), the 3D population inversion An 3D = {n 2 -n 1 ) can 
be found. Substituting this back into equation (III.D.34), the photon density is determined 
to be 


photon 


' photon 'Hi 

A\e\L 




where the internal quantum efficiency is 
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and the threshold current satisfying equation (III.D.34) is 
d|e|Lr.1 
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(III.D.35) 


(III.D.36) 


(III.D.37) 
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The power emitted by the laser can be found from P = ti (oV c n p hoton / r m , substituting for 

fiphotoni ^m-, and r 


P = Vi 01 m ^(/-4) (III.D.38) 

I e | 

The first three terms in equation (III.D.38) grouped together, is the differential quantum 
efficiency, which is the number of photons emitted per electron. The internal quantum 
efficiency (defined above at some biases) shows the output power is reduced by an 
imperfect (less than unity) injection efficiency rj and a finite lower lasing state lifetime 


T\. 


Below threshold, stimulated emission can be ignored because there are few 
photons within the cavity, Vn p hoton ~ 0. Solving for n 2 from the equation (III.D.32) and for 
m from the equation (III.D.33), a relation for the population inversion is found, noting 
that «2 = rj(J/\e\L)T 2 . Writing this in terms of the 2D population inversion 


A J 

Af hD=Vr-J. 




v 


l 21 J 


-0- ? ;)rrT 
e 


(III.D.39) 


61 



This equation allows for estimation of the threshold current. This simply shows it is 
necessary to have T\< r?i for a population inversion, that it is desirable to have a perfect 
injection efficiency 77 = I and a long upper state lifetime. Improving any of these 
parameters would serve to improve the population inversion. Because rate equations are 
not well suited for analyzing the electron transport in QC structures, in the next section, it 
will be discussed how Monte Carlo simulations can be used to model these structures. 
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IV. MONTE CARLO SIMULATION 


Since the lifetimes of the electronic states within QC structures are not constant, 
rate equations are not well suited for analyzing complete structures because they would 
be nonlinear. Thus, one approach that is more suited for this type of transport analysis is 
that based on Monte Carlo simulations. This approach does not rely on any assumptions 
about the electron distributions and can handle density and temperature dependent 
scattering mechanisms. By analyzing the complete structure, subband populations and 
electron distributions for all energy states can be detennined. Quantities such as current 
density (which is measurable) and other quantities such as gain can be found. This is 
important because while having a large oscillator strength is favorable, so is having a 
long upper state lifetime to maintain the 3D population inversion A« 3 d, as can be seen 
from the equation for peak gain = 2 e~h An^/m eva>FWHAd\(J\d/dx\i)\ (equation 
(III.D.27)). These are often competing characteristics. 

A. ELECTRON TRANSPORT IN QC STRUCTURES 

As has been discussed, the electron transport in entire QC structures is not suited 
for rate equation analysis, and the large number of scattering rates that must be computed 
would further complicate attempting to use such an approach. From the transport analysis 
of QC structures, we ultimately seek to find the populations and electron distributions of 
all states in the structure. The time evolution of the electron distribution functions is 
described by 

%=i -f v )~ (‘ - /«•)] (iv.a.i) 

at k ^' 

where S is the total scattering rate (sum of all mechanisms). The Monte Carlo method, 
applied to analyzing the electron transport in QC structures, consists of tracing a 
sufficient number of electrons throughout the structure over time, such that the results 
will be a good representation of the electron ensemble within the device. It relies on the 
use of random number generation, to appropriately choose from the probability 
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distributions encountered in this type of a simulation. There are several equivalent ways 
to formulate this technique. 95 ' 99 In what follows, it will be discussed how we chose to 
implement this technique. 


B. MONTE CARLO METHOD 

In general, the motion of a charged particle is due to the presence of fields, 
electric and magnetic, and scattering events. The time durations between scattering 
events, the scattering times, historically has been called the free flight times. During its 
time evolution, the electrons subband and momentum change, and so does the total 
scattering rate for all mechanisms. Consequently, the free flight times will change too. In 
order to be able to trace electrons throughout a structure, we must develop a method for 
stochastically generating these flight times. For example, consider the simple case of an 
electron moving under the influence of an applied electric field in the xi direction, the k\ 
component changes as k\=k\^o-\e\lh\E\\t, while the other momentum vector components 
remain unchanged. At the end of a free flight, the electrons momentum and energy are 
updated, and the electron is then scattering into its next state. 

In QC structures, an electron is characterized by its subband and components of 
the in-plane momentum vector k. Because the wavefunctions already contain the effects 
of the applied electric field, there is no need to include the momentum component that is 
in the growth direction. Next it will be discussed how to stochastically generate these free 
flight times for the special case when the scattering rate is constant. As will be shown, 
this fonnulation is still applicable to the more practical case where the scattering rate 
changes over the simulation time, as will be necessary for modeling the electrons within a 
QC structure. 

1. Flight Times 

Considering the scattering rate out of a subband to be constant T 0 = 1 It, an 
electrons probability to undergo a scattering event in a duration dt is T 0 dt, and dn = 
-T 0 ndt. This simple equation can be solved to find how the population changes with time. 
n(t) = n( 0)e~ r °‘ (IV.B.l) 
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The probability that the electron will scatter between t and t+dt, depends on the number 
of electrons in the subband and the scattering rate. This relation may be written as 
dtP(t) = dtr o e~ r °‘ (IV.B.2) 

Q7 

where P is a random number from 0 to 1. Using the direct continuous technique to 
stochastically generate the flight times, from equation (IV.B.2) 


r = 


J dtP(t ) 


(IV.B.3) 


o 

where r is a random number between 0 and 1. The solution for this equation is 
r = -e~ r ° ?/ +1 (IV.B.4) 

from which the flight time tf= -1/ T 0 ln(T-r) is solve for. Since r is a number between 0 
and 1, so is (1 -r), and it is equivalent to write 


t f = ~Y Xn ^ 


(IV.B.5) 


This expression allows us to stochastically generate the flight times for a constant total 
scattering rate. 

As discussed earlier, the scattering rates will not be constant. Nevertheless, this 
formulation is still applicable to the case where the total scattering rate changes over the 
simulation time. To deal with this, we will use a clever technique formulated in the given 
references, 95 ' 96 called self scattering. Self scattering is the difference between the constant 
scattering rate and the actual rate. The actual scattering rate is not constant and will 
depend on, for instance, the scattering mechanism, in-plane momentum vector k, subband 
of the initial electronic state, electron densities, and temperatures. 

r o = T seif +T(k) (IV.B.6) 


The reason for the inclusion of this self term, is that it simplifies the calculations that 
follow by allowing us to treat the scattering rate as being constant. The self term will 
eventually be removed, so that it will have no affect on the final results of the simulation. 
For this technique to work, T 0 must always be chosen greater than the maximum actual 
scattering rate T(k). The penalty paid, is that at the end of every free flight time, it must 
be checked to see if an actual scattering event or a fictitious self scattering event 
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occurred. If a fictitious self scattering event occurred, no real scattering event occurred 
and the electrons momentum is updated and a new flight time is generated. If on the other 
hand it is detennined that a real scattering event occurred, the electrons energy and 
momentum are updated and it is then scattered into its new state, where its energy and 
momentum are updated for this new state. This process of checking to see if a real 
scattering or self scattering event occurred, can generate a considerable amount of 
oversampling (self scattering events). Despite the oversampling that occurs due to the 
inclusion of the self tenn, it is much better than the alternative which would require 

QO 

changing equation (IV.B.2) to following integral equation (which is in the exponent). 

t 

dtP(t) = dtT(k,t)e 0 (IV.B.7) 

This equation is not suited for efficient repeated numerical calculations, and for this 
reason the self scattering technique is used instead. 

To use equation (IV.B.5) to generate the flight times, along with the self 
scattering technique, T 0 must be found for the simulation. There are different ways to 
implement this, we choose to use the maximum rates defined below, summing over all 
mechanisms (denoted by m) the maximum scattering rates T maXi m 

r„=Z r — (IV.B.8) 

m 

where 

r max, m = max , max ’ k [rW( k )l (IV.B.9) 

/ 

for our initial mapping scheme in the selection process that will follow. This will thereby 
introduce an additional self tenn in the selection process, 100 that will also have to be 
removed using a straight forward rejection technique. This provides a convenient way of 
implementing simulations with a large number of scattering rates, with the alternative 
method being to choose the maximum rate as a function of k value. This method of 
finding T 0 is illustrated graphically in Figure IV.B.l. 
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Figure IV.B.l. Graphical illustration for finding r o , for the Monte Carlo method. 

2. Selection Process 

Now that r 0 can be found and the flight times can be generated, we must proceed 
with the selection process that determines the scattering mechanism (or if a self scattering 
event occurred), the final state, and final energy and momentum of the followed electron. 
Single electron scattering mechanisms are all handled very similarly, and electron- 
electron scattering is a special case requiring more parameters to be determined. It will 
first be describe how the selection process works for the single electron scattering cases, 
and then for the electron-electron scattering case. Beginning with a new free flight, a 
random number r 0 is generated between 0 and T 0 and mapped onto the different 
maximum scattering mechanisms rates, using the following bin walling procedure, 97 ' 101 
where the mechanism m = x is selected if 

Z r — «.i Sr o (IV.B.10) 

m= 0 m= 0 

If no self scattering occurs and mechanism m = x is selected, the final subband state /= y 
is chosen by generating another random number r\ between 0 and r max m=x , and mapping 
as follows 

< r, (1V3.11) 

/= 0 /= 0 
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Once the mechanism and final state subband have been determined, the final momentum 
vector kf = kf{ cos6*xi+ sin#x 2 ) must be found. As an example, the magnitude kj is easily 
found from the conservation of energy Etj = Etj or for LO-phonon scattering Ej,, = 
E t j±E[X), where E T = E n +h k Urn. The angle 6 can be found by using another Monte 
Carlo technique. First the angle is randomly chosen from a uniform distribution between 
0 and 2n. The self tenn is removed using a rejection technique, which is an iterative 
scheme keeping the choice if 

r<^~ (IV.B.12) 

P 

max 


and otherwise rejecting it. P(9) is found from the corresponding scattering mechanism 
equations by 


P{0) = 


K 

\dOK 


(IV.B.13) 


where K is the kernel inside the scattering rate integral. This selection process and how it 
unfolds is graphically illustrated in Figure IV.B.2. 
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Figure IV.B.2. Graphical illustration of the Monte Carlo selection process. 

For electron-electron scattering, many more additional tenns must be found. The 
subband i, k, and 6, of the initial followed electron (three parameters) are known, and 
these three parameters for the second or parent electron to scatter with must be found, as 
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well as for the two final state electrons. In our implementation, this is done by bin walling 
in a similar fashion as discussed earlier, to find the ijfg process (a mapping look up 
scheme is used). A parent electron is then randomly chosen from the other electrons 
within subband j of the ensemble, which determines j, kj, and Oj. Since/and g are known, 
only four parameters remain to be found (Ay, Of, k g , and 0 g ). Only one of these really needs 
to be found though, as the others are dependant by the conservation of energy and 
momentum. We chose ^randomly from a uniform distribution between 0 and 2n, then 
using equations from section III.B (equations (III.B.69) et cetera), the other three 
parameters can be found. The self term can be removed by solving for q\\ from equations 
(III.B.52) or (III.B.71), and using the rejection technique of equation (IV.B.13). 

C. IMPLEMENTATION 

A flow chart describing the Monte Carlo simulation is shown in Figure IV.C.l. 
The algorithm consists of calculating the scattering rates and the main Monte Carlo 
selection process previously described. Because a QC structure is periodic and charge is 
conserved, when an electron with in-plane wavevector k is scattered into another section, 
an electron is introduced with the same k vector into an equivalent state in the section 
being modeled. Typically, the center section of a three section structure is modeled. 
Initially the electrons are distributed thermally at the lattice temperature according to 
Fenni-Dirac statistics, in one or more energy levels. In theory, the scattering rates need to 
be recalculated every time an electron scatters to a new state. However, in practical 
implementation of this algorithm, the scattering rates need only be calculated after an 
evolved time step duration At, and this duration of time is chosen such that the scattering 
rates do not change too much. For example, a At ~ 0.1 psec might be needed for quickly 
changed transient times. So in practice, each electron in the ensemble is traced for a time 
step duration At, then the scattering rates are recalculated based on the new subband 
populations and electron distributions. In this way, the scattering rates are consistently 
solved for using the correct subband distributions, as the time evolution of the ensemble 
unfolds. This allows for state blocking to be taken into account, and ensures that density 
and temperature dependant scattering rates are correctly modeled. The number of 
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particles must be chosen enough for accurate calculation, and typically ~ 10 4 is adequate 
for most simulations. Note that it is not that important how many electrons or the sheet 
density each particle represents, only that there are enough for accurately statistically 
modeling the ensemble within the structure. 


Ensemble 

loop 



Figure IV.C.l. A flow chart for the Monte Carlo simulation of QC structures. 
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Because the scattering rate calculations are the time consuming part of the 
simulation, we choose to calculate and store the rates before each time duration. That 
way, the rates are only calculated once for each time step duration. The Monte Carlo 
algorithm is continued until it is detennined that state populations have reached steady 
state and are not changing. This is done by making sure the variance of the populations is 
below some minimum acceptable value. Once steady state is reached, the electron 
densities and temperatures are output, and other parameters such as various scattering 
rate lifetimes, current density, and gain are easily found (all scattering events are stored). 
As an example of a typical simulation, shown in Figure IV.C.2 is a Monte Carlo 
simulation of QC structure D619F10E which has six energy levels within one section (the 
details of this structure will be discussed later), where all electrons were initially 
distributed in state 2 at the lattice temperature of F= 25 K. 



Figure IV.C.2. Monte Carlo simulation of structure D619F10E at a bias of 53.6 
mV/section, with all electrons initially distributed in state 2. Inset: Two sections of the 
structure are shown, with one section outlined. 

Under the transport picture described by this type of Monte Carlo simulation, 
where the scattering between electronic states is determined by the spatial wavefunctions 
solved via Schrodinger’s equation, coherence effects are not taken into account. For 
example the peak current through a thickened barrier is not affected, only its transmission 
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sharpness decreases, and this is not accurate. This implies that the scattering is not a fully 
coherent process and scattering between weakly coupled states with a small anticrossing 
is mostly a noncoherent process, due to dephasing scattering that effectively interrupts or 
scrambles the phase coherence. No phenomenological dephasing parameter is introduced 
in our simulations. This can sometimes lead to optimistic results as well as an 
overestimation of the scattering between states with a small anticrossing, and hence an 
overestimation of the current density when aligned with weakly coupled parasitic current 
channels. The scattering between more than one section is taken into account in step well 
structures because of the nature of the step well QC structure (which potentially can 
further lead to an overestimation of the scattering between weakly coupled states). The 
scattering mechanisms that are included are LO-phonon, electron-electron, impurity, and 
interface roughness scattering. These scattering rates are calculated via Fermi’s golden 
rule. All rates from these simulations are net rates and include backscattering. State 
blocking and screening are also included in the simulations. Next we discuss waveguide 
structures commonly used in QCLs for the confinement of emitted photons. 
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V. WAVEGUIDE AND RESONATOR 


QCL devices are fabricated by growing N repeated sections for the active region, 
and then processing the wafer sample into a waveguide to form the laser resonator. In 
these sections, waveguide parameters of interest and the two common types of waveguide 
structures, surface plasmon and metal-metal waveguides, are discussed. A third structure, 
a metal-patterned substrate waveguide was also investigated. 


A. QCL RESONATORS 

The waveguide and resonator are important, because the losses essentially set the 
threshold condition. At threshold, the gain equals the losses. The threshold gain can be 
found by considering a one round trip inside the resonator as follows 


I 0 e 2a ' L R 1i R l fi 2T *** , * L = I 0 


gmn th 
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--Un(f? ( ,i?i) + « H 

V 1L 
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(V.l) 


where I Q is the initial intensity, a m = -\/2Lki(R 0 Ri) is the mirror loss, a w is the waveguide 
loss, L is the length of the resonator, and R is the intensity reflectivity at a facet. The 
approximation that the mode gain = Txgain has been used. The confinement factor T, is 
defined as the ratio of the guided energy inside the active region to the total energy of the 
wave 

Rej Jd 2 xS-ii 

p Inactive region 

Rej J d 2 x S • n 

I all space 

where S = lAExH*. 

Since the waveguide mode is attenuated, the propagation constant is complex 
P= P'+iP", noting that a w = 2 P". In terms of the electric field and considering the case 
where the field facet reflectivity r = r m e'^ is the same for both facets, the following can 
be written 
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(V.3) 


In general, it is seen from the phase terms that the Fabry-Perot resonance condition is 
found P'L = njz-(j). By ignoring <f>, which is often done for optical systems, at resonance 
the usual condition is found L = nXJ2. Next, the types of waveguides used for QCLs will 
be discussed. 


B. SURFACE PLASMON AND METAL-METAL WAVEGUIDES 

Two types of optical confinement structures for THz QCLs have been used, the 
semi-insulating (SI) surface plasmon and the metal-metal waveguides. 1,32 102 ' 104 A third 
type of waveguide similar to the metal-metal waveguide, the metal-patterned substrate 
waveguide, was also analyzed (Figure V.B.l). 



(a) (b) (c) 


Figure V.B.l. QCL waveguides, the (a) surface plasmon, (b) metal-metal, and (c) metal- 
patterned substrate waveguides. 

The surface plasmon waveguide (Figure V.B.l(a)) consists of a metallic top and 
two side contacts, and a thin n+ doped semiconductor layer or plasma layer (where the 
real part of the pennittivity is negative) sandwiched between the active region and the SI 
GaAs substrate. Gaps are present between the sides of the active region and the metallic 
lateral contacts. This type of structure is easier to fabricate than the metal-metal 
waveguide since it does not require substrate removal and wafer bonding. Since the fields 
extend substantially below and outside the active region, the confinement factor is less 
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than that of the metal-metal waveguides. However, the advantage can be improved output 

13 

coupling for higher power applications. 

The metal-metal waveguide consists of two metallic layers above and below the 
active region (Figure V.B.l(b)) which is above a n+ GaAs substrate. The fields of the 
mode are usually well confined with often a near unity confinement factor, which can 
lead to lower threshold gains. The waveguide width may often be made smaller, than the 
surface plasmon waveguide counterpart, with acceptable performance. This can allow for 
reduction in the device area. The metal-metal waveguide is also better suited for higher 
operating temperatures. The disadvantage is that the output coupling is worse due to the 
high reflectivity of the subwavelength rectangular aperture, and the radiated beam pattern 
is much more fan shaped. At THz frequencies, we are in a waveguide region in which the 
Fresnel reflection coefficient is not valid to use for the reflectivity for the metal-metal 
waveguide. The high reflectivity is an advantage in the sense that a high reflectivity (HR) 
coating on one facet end is rarely needed with this type of waveguide. Although perhaps 
difficult to fabricate, others have shown improved output power in metal-metal 
waveguides by micromachining a matching horn antenna on the facet end. 105 Both the 
surface plasmon and metal-metal configurations use surface plasmons attached to the 
contact layers. 

The waveguide losses and reflectivities for the waveguides, are found using finite 
element method (FEM) solvers. ’ The Drude model is used to determine the 
material parameters, in which the complex permittivity can be found from 


/ \ • ^3D^ 

£(C0) = £ b +l - -j- -T" 

moo —ico 

U J 


(V.4) 


where \e\ is the charge of the electron, m is the effective mass of the electron, and co is the 
frequency. In our simulation model, the active region was nominally taken as GaAs and 
Sb was set to 12.4 corresponding to the f=0K value. The relaxation times r, were set to 
0.1 ps for the n+ layers, 0.5 ps for the active region, and 0.06 ps for the metallic gold 
layers. 108 ' 110 For the active region a background carrier concentration of 2xl0 1:i cnT 3 was 
assumed, and a concentration of 5><10 cm was used for the metallic gold layers. The 
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plasma and contact layers concentrations vary by design, and are listed in those particular 
sections. From Ampere’s law V *H = J+dD/dt = gE-/7uD and equation (V.4), it is seen 
the complex permittivity can alternatively be entered as a real valued permittivity Sb and a 
complex conductivity a = n^e T/m(\-io)z). Additionally, it can also be entered as a real 
valued effective permittivity s' and a real valued effective conductivity ox” since V xH = 

ox "E -ico sE, noting that s' can be a positive or negative value. By considering 0 ) in the 

2 1/2 

high frequency limit, the plasma frequency may be defined as co p = {nyae~lz\,m ) . Since 
ck= (o) - 0 ) p ~ ) it is seen that when ox co p , k is purely imaginary and waves incident on 
the plasma are reflected and waves inside fall of exponentially with distance, thus the 
medium behaves more like a “metal.” For a>> (o p transmission occurs and the medium 
has more of a lossy “dielectric” characteristic. 

A waveguide similar to the metal-metal waveguide, the metal-patterned substrate 
waveguide, was also briefly investigated (Figure V.B.l(c)). The idea was if MBE growth 
of the active region could be accomplished on top of a patterned substrate, this could 
avoid having to wafer bond while at the same time keep a structure capable of 
maintaining a high confinement factor. It was ultimately determined that the defects 
introduced by the growth process over the patterned substrate, would likely be 
unacceptable for high quality MBE growth. Nevertheless, the results of the analysis will 
be briefly mentioned, as similar type of 3D analysis could be useful for mode suppression 
waveguides. In such a waveguide, a metallic layer would be deposited on the top surface 
of the active region, and the pattern would also have to be dense enough to keep the 
modes confined while also not substantially increasing the waveguide loss. Our analysis 
showed for the guides we modeled, compared to metal-metal waveguides, that the 
confinement factor was affected little by square patterns from 1 by 5 jum (1 fjm thick 
strips with 5 /an spacing between the strips) to 2 by 2 //m, but the waveguide loss was 
increased for the less dense patterns. It was also shown for a 2 by 2 fjm pattern 
waveguide, that the waveguide loss could be made lower using a n+ substrate (as 
opposed to a SI substrate with a thin n+ layer) for frequencies beyond about 2.5 to 5 THz. 
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VI. THz LO-PHONON QC STRUCTURES 


Two different LO-phonon approaches were investigated, the new step well QC 
structures we proposed, 40 ’ 111 ' 113 and a conventional square well design. Approximate rate 
equation analysis was initially perfonned on the step well structures, and the structures 
were analyzed using Monte Carlo simulations. The analysis indicated the injector of one 
of the step well structures could further be improved, and simulations on a modified 
structure showed a more uniform gain over a range of biases. Simulations were also 
performed on a high power square well LO-phonon design. Our analysis was compared 
to experimental measurements, and found to be in reasonable agreement. 

A. STEP WELL QC STRUCTURES 

Even though it is possible to arrange the energy levels in a step well such that the 
radiative transition and LO-phonon transition are within the same step well, this does not 
mean it will necessarily be useful for a laser application. For the device to lase, a 
population inversion must be kept for sufficient gain to overcome the losses. While 
intrawell radiative transitions can have large overlap of the electron wavefunctions that 
yield large oscillator strengths, there can be a trade-off between the oscillator strength 
and upper state lifetime, as the scattering from the upper radiative state to the lower state 
is usually increased with increasing oscillator strength in three-level LO-phonon 
structures. This increase in oscillator strength is one of the reasons for studying these step 
well structures, as increasing the oscillator strength can increase the gain and hence the 
output power of a device. As will be shown, the section length will remain about the 
same as other LO-phonon designs. 

Initially, an approximate rate equation analysis was perfonned on the step well 
structure shown previously in Figure II.B.l(b). 40111 The scattering rates (W 2 \ and W \ 0 ) 
must be calculated for the step quantum well structures in order to estimate whether a 
population inversion is likely between states 2 and 1. The simplest rate equation analysis 
allows us to write the population of the middle state 1 (n\) as 
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- L = n 2 W 2l -n l W l0 (VI.A.1) 

at 

where n 2 is the population of state 2. The necessary and intuitive condition at steady state 
for a population inversion to exist ( n 2 > «i) is JVio > W 2 \ (in terms of corresponding 
lifetimes rio < t 2 \). In general, the scattering rates are a combination of all possible 
scattering mechanisms, i.e., electron-phonon, electron-electron, impurity, and interface 
roughness scattering. In this approximate rate equation analysis, we estimate the 
scattering rates by taking into account the LO-phonon and electron-electron scattering 
rates. As will be seen later from our Monte Carlo analysis, the electron temperatures can 
be ~ 50 to 100 K or higher than the lattice temperature. For the rate calculations that 
follow, the lattice temperature will be taken to be 25 K and the electron temperature to be 
100 K. 

Figure VI.A.l shows the mean LO-phonon scattering rates for the 4.3 THz step 
quantum well previously shown in Figure II.B.l(b) (and shown again inset), calculated 
for initial state populations from lxlO 9 to lxlO 10 cm -2 for T/ allice = 25 K and T e / ectmn = 100 
K. All of the scattering rates are essentially constant and independent across the 
populations computed. As expected, the 2 to 1 transition has the smallest mean scattering 
rate and hence the longest lifetime, because the energy spacing of 17.9 meV between 
states 2 and 1 is below the LO-phonon energy spacing of about 36 meV. Only the 
electrons with sufficiently high in-plane kinetic energy (~ 36-17.9 meV = 18.1 meV or 
higher) can participate in the scattering process, while others are forbidden. The 
maximum scattering rate, which occurs when the electron has just enough in-plane 
energy such that the energy difference is exactly at the resonant LO-phonon energy Hcolo, 
is also shown in Figure VI.A.l. It can be seen that the mean scattering rates for the 1 to 0 
transition are faster than the 2 to 1 transition, and are even faster than the 2 to 1 
maximum rate. Also included is the parasitic mean scattering rate for the 2 to 0 transition. 
This parasitic lifetime is 1.8 picoseconds (psec), which is relatively short. These 
scattering rates show that the mean lifetimes are such that t 2 \ = 5.6 psec > rio = 0.37 
psec. This indicates that the LO-phonon rates are favorable for keeping a population 
inversion. Of course, the population inversion is affected by all scattering mechanisms, 


78 



not just by the electron-phonon scattering rates. Next we will calculate the electron- 
electron scattering rates and see how incorporating these rates may affect the population 
inversion analysis. 



2 4 6 8 10 

n(10 9 cm' 2 ) 


Figure VI.A.l. Mean LO-phonon scattering rates for the 4.3 TFIz step quantum well 
(shown inset) as a function of initial state populations for Ti attice = 25 K and T e i ectron = 100 
K (solid lines). For the 2 to 1 transition the maximum scattering rate (which is not an 
averaged rate) is also shown (dashed line). The lifetimes are listed (reciprocals of the 
essentially constant scattering rates). 

The electron-electron scattering rates were calculated for the same 4.3 THz 
structure as for the LO-phonon analysis using T e i ectron = 100 K for initial state populations 
of lxlO 9 to lxlO 10 cm 2 . Figure VI.A.2 shows the electron-electron scattering rates for 
the 2 to 1 and 1 to 0 transitions that pertain to the simplified rate equation analysis 
already discussed. It is seen that the dominant scattering rate is the 2010 process. Unlike 
the earlier LO-phonon analysis, the electron-electron scattering rates are not constant but 
rather as expected, increases with carrier concentration. The next highest scattering rates 
are for the 2221 process, and the lowest scattering rates are for the 1110 process. The 
total electron-electron scattering rates are the sum of the above processes, taking into 
account the number of electrons that change energy levels. Figure VI.A.3(a) shows the 
total electron-electron scattering rates for the 2 to 1 and the 1 to 0 transitions, based on 
the processes for the simplified rate equation analysis. It is clear that the 2 to 1 scattering 
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rates are higher than the 1 to 0 scattering rates for all the initial carrier populations. Thus 
the lifetimes are such that the electron-electron scattering process by itself is not 
favorable for the purpose of keeping a population inversion. 

Although our analysis so far has been in the framework of the simplified rate 
equation analysis and has considered only processes within that framework, it is worth 
mentioning that in general a system with N number of energy levels will have 4 N electro¬ 
electron scattering processes. Some of these processes do not directly by themselves 
affect subband populations, i.e., intrasubband iiii scattering rates where no net change 
occurs in the number of electrons in the subband. To see what the scattering rates are for 
some of the additional processes, computed in Figure VI.A.3(b) are the mean electron- 
electron scattering rates for the additional 2110, 2210, and 2220 scattering processes. It is 
seen that of these three, the 2110 process is dominant and involves a net change of losing 
one electron from state 2 and gaining one electron in state 0, a process by which 2 —> 1 
and 1 —> 0. The scattering rate due to this process is around an order of magnitude less 
than those considered earlier. 



Figure VI.A.2. Mean electron-electron scattering rates for the 4.3 THz step quantum well 
as a function of initial state populations for T e i ectron = 100 K. The various scattering 
processes for the 2 to 1 and the 1 to 0 transitions for the simplified rate equation analysis 
are shown. 
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Figure VI.A.3. Scattering rates for the 4.3 THz step quantum well for T e i ectron = 100 K. (a) 
The total electron-electron scattering rates for the 2 to 1 and the 1 to 0 transitions, based 
on the processes for the simplified rate equation analysis, as a function of initial state 
populations, (b) The mean electron-electron scattering rates as a function of initial state 
populations for T e i ectron = 100 K, for some additional scattering processes for the 2 to 1 
transition that are not within the framework of the simplified rate equation analysis. 


The analysis indicates that the electron-phonon scattering rates are around an 
order of magnitude larger than the electron-electron scattering rates. The scattering rate 
analysis of the biased step quantum well indicates that this type of structure is likely 
capable of keeping a population inversion for the temperatures and range of carrier 
concentrations considered. 

In order to be used in a QC structure, the step quantum well must have an injector 
section to feed the next section of the cascade. One potential implementation is illustrated 
in Figure VI.A.4 where one section is outlined. 40 The double barrier injector section was 
designed to maintain the upper state lifetime, while allowing the electrons in the lower 
states to resonantly tunnel to the upper state of the next adjacent section. The conduction 
band profile was solved by using a self consistent solution to Schrodinger’s and Poisson’s 
equations. The step well has been shaped to spatially separate the upper state and lower 
states as much as possible, in an effort to reduce scattering between those states. The 
radiative transition takes place between states 4 and 3, where E 43 = 16.6 meV (4 THz or A 
~ 75 //m). The LO-phonon assisted transition takes place from state 3 to the lower triplet 
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states (2, 1, 0), since £31 = 36.6 meV which is near hco L o. Scattering between the triplet 
states ( 2 , 1 , 0 ) and injection into the next adjacent sections upper state is to take place. 
The step well has been arranged so that state 4 is above the step and state 3 below the 
highest point, in an effort to spatially separate the wavefunctions for efficient injection to 
the upper state 4 and for reducing the parasitic injection to state 3. Thus, the step barrier 
helps to reduce unwanted injection into the middle state, which is a problem with other 
square well LO-phonon structures. 



Figure VI.A.4. Conduction band profde of a proposed step well QC structure found using 
a self consistent solution to Schrodinger’s and Poisson’s equations. The step well and 
resonant tunneling double barriers of one section are outlined. Beginning with the left 
injector, the Al x Gai x As layers compositions are 0.143/0.035/0/0.143/0/0.143/0 and 
thicknesses in mn are 4.3/20.9/13.5/1.7/9.6/2.5/7.6, with £43 = 16.6 meV (4 THz) and £31 
= 36.6 meV. The applied field is 9.5 kV/cm. The 9.6 mn well is doped 2.9x 10 16 cnT 3 , 
which corresponds to a sheet density of 2 . 8 X 10 10 cm -2 . 

The mean electron-phonon scattering lifetimes for this structure are 
* 3 ->( 2 ,i,o) = 0-4 psec, for the upper to lower states scattering Z 4 _>( 2 ,i,o) = 1.4 psec, and 
r 43 = 7.9 psec. Again, because these are intrawell transitions, the scattering rates 
computed are relatively fast. The radiative transition is also an intrawell transition, and 
the overlap of the wavefunctions is reasonably good, which gives an oscillator strength of 
about 1.34. This oscillator strength is higher than the highest reported LO-phonon 
designs (~ 40% increase), and the section length is about the same, i.e., about half that 
of typical bound to continuum designs. 
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To get an approximate estimate of what might be expected for this type of 
structure, we can estimate from the gain equation (equation (III.D.27)) the population 
inversion. If we assume a total loss of ~ 40 cm 1 (which is reasonable for surface plasmon 
waveguides) and assume a FWHM ~ 4 meV, the estimated 2D population at threshold is 
~ 2.5><10 cm . This is within the carrier densities used in our analysis. The current 
density can be estimated from A« 2 d ~ (J/e) 1 - zV243 ) and is calculated to be J ~ 760 

A/cm 2 . 

It was concerned that the triplet of lower states spacing was too close to the 
radiative energy spacing and could provide a means for unwanted absorption. Further, 
more energy separation was wanted between the lower states and an unused state in the 
next adjacent section. For this reason, it was decided to focus on single well injector 
designs. There were a number of different step well structures that were considered, some 
of which are shown in Figure VI.A.5. Other structures (not shown), that resonantly 
tunneled the middle and lower states, step well designs without all three energy levels in 
the same well, and some structures that would also have required growth with more than 
two A1 compositions in an Al x Gai_ x As heterostructure, were considered as well. Due to 
the difficulty in growing a sample of repeated sections with more than two A1 
compositions in most MBE chambers, structures requiring only two A1 compositions 
were favored. 





Figure VI.A.5. Some of the step well QC structures considered, (a) one step, one well 
injector, (b) two step, one well injector, and (c) one step, two well injector. 
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By using a single well injector, the anticrossing between the lower ground states 
can be kept much smaller than the optical transition energy spacing (as there will only be 
a doublet of states), thereby reducing the possibility of absorption. Initially a single step 
well and one well injector (Figure VI.A.5(a)) structure was analyzed. To avoid a potential 
parasitic current channel (from alignment with an unwanted energy state in the next 
adjacent section), the injector well was narrowed so that the only energy state from the 
well was that of the lowest ground state. This increased the sub-picosecond LO-phonon 
lifetime of the middle state to ~ 1 psec which was too slow. Rather use a wider injector 
well to keep the sub-picosecond LO-phonon scattering lifetime and also introduce a 
second state from the well which could potentially produce a parasitic current channel, it 
was chosen to utilize a two step well and one well injector as shown in Figure VI.A.5(b) 
(where the two steps are of equal A1 composition). This structure, QC structure 
D619F10D, will be discussed in the next subsection. The disadvantage, however, is that it 
introduces two additional growth layers (7 growth layers as compared to only 5 growth 
layers for a single step well one well injector design). Nevertheless, even with a two 
equal A1 composition step well, this structure is quite simple having only two wells. 
Albeit the first being a step well. It also has only one additional layer compared to the 
simplest square well LO-phonon QC structure/ Note that in principle, a step well QC 
structure does not need an injector well (in which it would be a one well design). 

1. Structure D619F10D 

The conduction band profile for QC structure D619F10D is shown in Figure 
VI.A.6. Electrons are injected into the upper state 3 and the radiative transition occurs 
between state 3 and middle state 2, where E 32 =15.1 meV (~ 3.6 THz or 2 ~ 82 /an). The 
LO-phonon assisted transition takes place from the middle state 2 to the lower doublet 
states (1, 0), since E 2 \ = 37.2 meV which is near ha>w- This ensures fast depopulation of 
state 2 via LO-phonon scattering, with a mean scattering rate o) = 0.68 psec (0.44 
psec band edge). The anticrossing between states 2 and 5' is relatively small due to the 
thick collector and injector barriers, and therefore should not be a primary path of 
scattering out of state 2. Though any electrons that transition from state 2 to 5' would 
serve to help the population inversion, and state 5' is near the resonant LO-phonon energy 
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spacing to upper state 3'. Injection from the lower doublet states (1, 0) is to take place 
into the next adjacent sections upper state 3’, and the injector anticrossing Eoy - 1 meV. 
The step well has been arranged so that state 3 is above the step and the state 2 below the 
highest point so as to provide spatial separation at the injector barrier for improved 
injection. 



Figure VI.A. 6 . Conduction band profde of structure D619F10D found using a self 
consistent solution to Schrodinger’s and Poisson’s equations. The step well and resonant 
tunneling well of one section are outlined. Beginning with the left injector, the 
Al x Gai_ x As layers compositions are 0.16/0.05/0/0.05/0/0.16/0 and thicknesses in nm are 
4.8/27.9/1.8/2/6.3/4.1/6.8, with f /32 = 15.1 meV (3.6 THz) and E 21 = 37.2 meV. The 
applied field is 10.5 kV/cm. The center 2 nm of the 6.8 nm well is doped to a sheet 
density of3.4><10 10 cm~ 2 . 

Since the transitions are intrawell in nature, the LO-phonon scattering is relatively fast. 
The radiative transition is also intrawell, and the overlap of the wavefunctions is 
reasonably good, giving an oscillator strength of - 1.03 at 9.9 kV/cm (53.2 mV/section). 

In order to determine the threshold gain required for lasing, two different 
commonly used waveguide structures, the surface plasmon and metal-metal waveguides 
were considered. The threshold gain was determined for both surface plasmon and metal- 
metal waveguide resonator configurations at the operating frequency. The analysis was 
performed assuming an active region thickness of - 10 //m (which corresponds to N ~ 
185 sections) and a guide width of 200 jum. This is because the threshold gain is worse 
(higher) for smaller width ridges. The threshold gain = (a m +a w )IY, where a m is the 
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mirror loss and a w is the waveguide loss. The waveguide losses and reflectivities for the 
metal-metal waveguide were found using finite element method (FEM) solvers. 108,109 The 
Drude model was used to determine the material parameters, and a background carrier 
concentration of 2><10 15 cm" 1 was assumed. Due to the difficulty in accurately modeling 
the waveguide with N thinly doped active region section layers, the doping in these thin 
layers was neglected. Thus, the actual waveguide losses may be higher. 

Considering the surface plasmon waveguide first, the side contacts spacing was 
set at 50 nm to keep the mode from coupling to avoid higher waveguide loss. The top 
contact plasma layer was 60 nm thick and doped at 5><10 18 cnT 3 . To determine a suitable 
lower plasma layer thickness, the threshold gain was computed as a function of lower 

IQ 1 o _ ■j 

plasma layer thickness for doping concentrations of 1><10 to3 x 10 cm . The mirror 
losses have been computed for resonator lengths of 1 and 2.5 mm with one facet of the 
waveguide assumed to have a FIR coating. Figure VI.A.7 shows that for the lowest 
threshold gain over the range of mirror losses considered (resonator lengths of 1 and 2.5 

1 O 1 O _T 

mm), the optimal lower plasma layer is one doped ~ 2x10 to 3 X 10 cm at a thickness 
of 400 to 500 nm. 



Plasma Layer Thickness (nm) 


Figure VI.A.7. Threshold gain for a 200 jun i ridge, 10 jum thick surface plasmon 
waveguide, as a function of lower plasma layer thickness for different doping 

_o # _1 

concentrations (cm ) and mirror losses (cm ). 
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The cross-section 2D mode intensity and ID profile mode intensity are shown in 
Figure VI.A.8 for a lower plasma thicknesses of t v i asm a = 500 nm, doped at 3x 10 18 cnT 3 . 
The confinement factor is T = 0.28 and the waveguide loss is a w = 3.3 cnT 1 for this 
waveguide. 
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Figure VI.A.8. 200 jui n ridge, 10 /an thick surface plasmon waveguide with t p i asm a = 500 

1 O _T 

nm doped at 3x10 cm , (a) cross-section 2D mode intensity, (b) ID profile mode 
intensity. 

Next, the threshold gain was computed for a metal-metal waveguide, where the 
top and lower contact plasma layers were 60 nm thick (10 and 50 nm layers doped at 

IQ 1 O _T 

5 x 10 and5 x 10 cm , with the lower contact also having a 3.5 nm LTG GaAs layer). 
For the metal-metal waveguides, we assume both facets are uncoated. We find the 
threshold gain to be 7.7 and 5.7 cnT 1 for mirror losses corresponding to resonator lengths 
of 1 and 2.5 mm respectively. The attenuation was calculated to be a w = 4.3 cnT 1 with a 
confinement factor T ~ 1. This near unity confinement factor is seen as evident from the 
cross-section 2D mode intensity graphed in Figure VI.A.9. Optimal guide parameters for 
a surface plasmon waveguide depend on having a relatively thick lower contact layer 
compared to the upper contact layer, while for a metal-metal guide both contact layers 
should be relatively thin to minimize waveguide loss. Because an etch-stop layer can be 
used for processing both surface plasmon and metal-metal waveguides from the same 
sample, a sample that has an optimal lower plasma layer thickness for a surface plasmon 
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waveguide will not be optimal for a metal-metal waveguide. As an example, the 
threshold gain was computed for a metal-metal waveguide with a top contact layer 60 nm 

1 o _o 

thick doped at 5x 10 cm , and for lower plasma thicknesses of t p i asm a = 150 and 500 nm 

1 o _o 

doped at 3x10 cm (and found to change very little for guide widths from 100 to 200 
/an). The threshold gain was calculated to be only 8.5 cnT 1 for t v i asm a = 150 nm and 
increased to 18.2 cnT 1 for t v i asm a = 500 nm (both with a resonator length of 1 mm). The 
thinner plasma layer thickness of 150 nm can be seen to result in an improvement for the 
metal-metal waveguide. Although the threshold gain is higher for a surface plasmon 
waveguide, the output coupling factor aj(a m +a w ) can be smaller for some surface 
plasmon waveguides which is better for high power applications, i.e., since the output 
power = VirjiaJtam+a^Nhco l\e\(I-I,h) from a single facet of an uncoated device. 



Figure VI.A.9. 200 /an ridge, 10 /an thick metal-metal waveguide, cross-section 2D 
mode intensity. 


a. MBE Growth and Fabrication 

The sample was grown by MBE using IQE foundry service, with N = 185 
sections, where the top and lower contact plasma layers were 60 nm thick (10 and 50 nm 

IQ 1 o _T 

layers doped at 5 X 10 and 5 x 10 cm , with the to be lower contact also having a 3.5 
nm LTG GaAs layer). The sample was processed by the National Research Council 
(NRC) (Canada) into metal-metal waveguides of 100 and 200 /an wide, and 1 mm 
resonator lengths, with both facets uncoated (Figure VI.A. 10). The fabricated devices 
showed an unusually high impedance (> 1 kO) at room temperature, compared to typical 
QCLs of about 10 Q. This likely indicates that there was either a growth problem, such as 
a high p-type background concentration during MBE growth (which depleted the electron 
concentration in the active region), or non-ohmic contacts generated during the 
processing. Initial electroluminescence measurements were made at T = 10 K (using a 
CTI Cryogenics Model 22 closed cycle He cryostat and a Nexus 870 FT-IR ESP Nicolet 
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spectrometer), and no light out was detected. Measurements on the sample showed little 
current through the device near the injector anticrossing (at -10.6 V, - 300 mA, which 
corresponds to 150 A/cm ), as illustrated in Figure VI.A. 11 for a 200 /an wide device. 
The I-V characteristics were measured on both the 100 and 200 //m wide samples out to 
about - 5 A, without the device going into NDR. Because the contact resistance was so 
high and no light was detected, there was either a growth or contact processing error with 
this sample. 




(b) 


Figure VI.A.10. Sample D619F10D. (a) Top view of the metal-metal waveguides, 
showing bond wires to the 100 and 200 /an wide by 1 mm long waveguides, (b) Cross- 
section of one of the 100 /an wide waveguides. 
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We are currently investigating the cause of the high impedance to determine the origin 
which will help us refine the growth and fabrication processes. Tests are being 
performed, using secondary ion mass spectrometry (SIMS) at NRC, to determine if 
correct doping concentrations were maintained in the sample growth. No conclusion on 
the feasibility of this device can be drawn from these measurements, as clearly there was 
a fabrication problem. 



Figure VI.A. 11. The measured I-V curve of sample D619F10D, for a 200 //m wide by 1 
mm long metal-metal waveguide device. 

b. Monte Carlo Simulation 

A Monte Carlo simulation of QC structure D619F10D was performed, to 
model the electron ensemble. As previously discussed, the scattering between more than 
one section is taken into account because of the nature of the step well QC structure 
(which potentially can further lead to an overestimation of the scattering between weakly 
coupled states). The scattering mechanisms that were included are LO-phonon, electron- 
electron, impurity, and interface roughness scattering. The scattering rates were 
calculated via Fermi’s golden rule. All rates from these simulations are net rates and 
include backscattering. State blocking and screening are also included in the simulations. 
The Monte Carlo simulations of the step well QC structure were perfonned with 7) atece = 
25 K. Shown in Figure VI.A. 12 is the calculated current density and gain (assuming a 4.1 
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meV FWHM) over the high gain bias region. The calculated gain is higher over much of 
this region than the threshold gains computed earlier for the waveguide resonator 
configurations analyzed. However, just beyond 50 mV/section a dip in the gain curve is 
observed. It was felt that a more uniform gain could be realized, by modifying the 
structure slightly to improve the injection efficiency. 




Figure VI.A. 12. Results of the Monte Carlo simulation of structure D619F10D, with and 
without including interface roughness scattering, (a) Current density as a function of 
applied electric field bias, (b) Gain as a function of bias. 

2. Structure D619F10E 

The dip in the gain curve seen from the Monte Carlo simulation of QC structure 
D619F10D, was likely a result of inefficient scattering into the upper state at those biases. 
In order to correct for this, a second structure was designed and analyzed, structure 
D619F10E where the injector well was thinned slightly and the collector well thickness 
increased slightly. With these changes, it was anticipated a more uniform gain over these 
biases could be achieved. The conduction band profile for the step well QC structure is 
shown in Figure VI.A.13. 113 Electrons are injected into the upper state 3 and the radiative 
transition occurs between state 3 and middle state 2, where £32 =15.2 meV (~ 3.7 THz or 
X ~ 82 //m). The LO-phonon assisted transition takes place from the middle state 2 to the 
lower doublet states (1, 0), since £21 = 37.9 meV which is near hcoio (~ 36 meV in 
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GaAs). This ensures fast depopulation of state 2 via LO-phonon scattering, with a 
scattering rate lifetime ^-qi.o) ~ 0.5 psec (band edge) for this structure. Again, the 
anticrossing between states 2 and 5' is relatively small (less than ~ 0.5 meV) due to the 
thick collector and injector barriers, and therefore should not be a primary path of 
scattering out of state 2. Though any electrons that transition from state 2 to 5' would 
serve to help the population inversion, and state 5' is near the resonant LO-phonon energy 
spacing to upper state 3'. Injection from the lower doublet states (1, 0) is to take place 
into the next adjacent sections upper state 3', and the injector anticrossing has been 
slightly increased to E$y ~1.2 meV. The step well has been arranged so that state 3 is 
above the step and the state 2 below the highest point so as to provide spatial separation 
at the injector barrier for improved injection. 



Figure VI.A. 13. Conduction band profile of structure D619F10E found using a self 
consistent solution to Schrodinger’s and Poisson’s equations, shown at an applied electric 
field of 10.4 kV/cm (56.3 mV/section). The step well and resonant tunneling well of one 
section are outlined. Beginning with the left injector, the Al x Gai x As layers compositions 
are 0.16/0.05/0/0.05/0/0.16/0 and thicknesses in mn are 4.6/27.9/1.8/2/6.3/4.7/6.8, with 
E 32 = 15.2 meV (~ 3.7 THz) and Ej\ = 37.9 meV. The center 2 mn of the 6.8 mn well is 
doped to a sheet density of 3.4x 10 10 cm -2 . 

Since the transitions are intrawell in nature, the LO-phonon scattering is relatively fast. 
The radiative transition is also intrawell, and the overlap of the wavefunctions is 
reasonably good, giving an oscillator strength of ~ 0.94 at 9.9 kV/cm (53.6 mV/section). 

The Monte Carlo method was used to model the electron ensemble of QC 
structure D619F10E, in the same manner as for the previous structure analyzed. 113 The 
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Monte Carlo simulations of the step well QC structure in Figure VI.A.13 were performed 
with Tlattice = 25 K. Table VI.A.I lists the subband populations and electron temperatures 
found at an applied electric field of 53.6 mV/section. 


n 

« 2 D (xlO 9 cm 2 ) 

Telectron (K) 

0 

16.1 

131 

1 

10.1 

132 

2 

1.53 

155 

3 

6.22 

127 


Table VI.A.I. Populations of the subbands ( E n ) and electron temperatures found from the 
Monte Carlo simulation of structure D619F10E, at an applied electric field of 53.6 
mV/section, with Ti attice = 25 K. 

The current density, gain, and population densities and the electron temperatures for the 
upper (£ 3 ) and middle (Ei) states are shown in Figure VI.A.14 over a range of applied 
electric field biases. The results are shown with and without including interface 
roughness scattering. Others have shown reasonable agreement with experimental 
findings by not including interface roughness scattering in their simulations , 114 and that 
agrees with our simulations of other square well structures. In these simulations interface 
roughness scattering was found to increase somewhat the current density and electron 
temperatures, and also slightly decreases the medium gain. The current density as a 
function of applied bias follows about what is expected, with the exception of the points 
between ~ 45 to 50 mV/section. This is where some scattering to parasitic current 
channels takes place. These high values are likely due to an overestimation of the 
scattering to the weakly coupled states. A small bump in the curve at ~ 38.4 mV/section 
is observed, where this is near the anticrossing between states 0 and 2'. The step serves to 
also minimize scattering between states 0 and 2 ’ at this intennediate bias. 
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Figure VI.A.14. Results of the Monte Carlo simulation of structure D619F10E, with and 
without including interface roughness scattering, (a) Current density as a function of 
applied electric field bias, (b) Gain as a function of bias, (c) The population densities and 
(d) electron temperatures for the upper (C 3 ) and middle (E 2 ) subbands as a function of 
bias. 


A peak gain of 86.9 cnT 1 is found (assuming a FWHM =4.1 meV) and the gain is 
~ 70 cm 1 or higher across much of the high gain bias region. This is due to the vertical 
transitions of the optical and LO-phonon transitions, as well as the high injection 
efficiency from the step. The vertical optical transition allows for a relatively high 
oscillator strength and fast LO-phonon scattering helps keep the middle state 2 depleted. 
Due to the step, injection loss to the middle state 2' is relatively small. The injection 
efficiency is found to be ~ 89%, defined here to be from the lower doublet states (0, 1) 
to the upper state 3'. The lifetime of the middle state due to LO-phonon scattering is 
found to be T 2 -^i,o) = 0.88 psec, with the ratio of r 3 ^(i,o)/z 2 -<i,o) from LO-phonon 
scattering varying from about 2.1 to greater than 3 across biases in the high gain region. 
The calculated gain is seen to be higher than the calculated threshold gains for both the 
surface plasmon and metal-metal waveguides analyzed earlier. 


94 


































B. SQUARE WELL LO-PHONON QC STRUCTURE 

As previously discussed, only square quantum well structures have been used in 
QCL designs. The highest power QCL, reported in the given reference, 10 utilized a four 
square well LO-phonon structure and operated at ~ 4.4 THz. This device represents a 
relatively optimized design in terms of the structure, wavelength of operation, and 
waveguide, to achieve a peak power of ~ 250 mW in a ~ 200 /an wide by ~ 1.2 mm long 
surface plasmon waveguide configuration. There were N = 183 sections grown, and the 

1 O _-5 

lower plasma layer was 400 mn doped at 3><10 cm . As with all conventional square 
well LO-phonon designs, it contained the usual first three square wells, and a fourth 
injector well similar to the first LO-phonon THz QCL. The third well is used to 
resonantly tunnel the lower lasing state and also arrange for a lower ground state. 
Because of this, as with other LO-phonon square well designs, there is a doublet of states 
at the lower lasing state. This structure will next be described in detail, and we will model 
it using our Monte Carlo code to compare the simulation results to the experimentally 
measured results. 

1. 4.4 THz Structure 

The conduction band profile for the 4.4 THz QC structure is shown in Figure 
VI.B.l. Electrons are injected into the upper state 4 and the radiative transition occurs 
between state 4 and state 3. The design was intended to emit near 4.7 THz, but was 
experimentally found closer to ~ 4.4 THz (/_ ~ 68 //in). The LO-phonon assisted 
transition primarily takes place from state 2 (as well as from state 3) to the lower doublet 
states (1,0), since £21 is near ha>w (~ 36 meV in GaAs). States 3 and 2 are intended to be 
coupled for fast scattering between those states. This ensures fast depopulation of the 
lower lasing state via LO-phonon scattering. Injection from the lower doublet states (1,0) 
is to take place into the next adjacent sections upper state 4'. This shorter wavelength of 
operation and hence larger radiative energy separation of the upper and lower lasing 
states, along with good coupling at the injector, aids in improved injection. 
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Figure VI.B.l. Conduction band profile of the 4.4 THz structure, shown at an applied 
electric field of 11.5 kV/cm (62.4 mV/section), with one section outlined. Beginning with 
the left injector, the Al x Ga[ x As layers compositions are 0.15/0/0.15/0/0.15/0/0.15/0 and 
thicknesses in nm are 4.8/8.2/1.7/6.8/4.0/16.4/3.4/9.0. The 16.4 nm well is doped to a 
sheet density of 3.1 X 10 10 cm -2 . 

The Monte Carlo simulations of the 4.4 THz QC structure were performed with 
Tiattice = 25 K. The current density and gain are shown in Figure VI.B.2 over the high gain 
bias region of the device. The results are shown with and without including interface 
roughness scattering, for roughness of a/2 (same as previously used for the step well 
structures analysis) and a/4. In these simulations, reasonable agreement is shown with the 
interpolated measured values 10 for the current density at most biases, with the largest 
difference being an overestimation occurring ~ 68 mV/section which is just beyond the 
injector anticrossing and near the anticrossing. This overestimation is likely due to 
not including coherent effects of dephasing. As with the previous step well structure 
analyzed, the effects of interface roughness scattering were found to increase somewhat 
the current density, and also slightly decrease the medium gain. It is seen that assuming a 
roughness of a/2 overestimates the current density, and that using a roughness of a/4 is 
almost negligible compared to simulations without including interface roughness 
scattering. Although one monolayer (ML) is sometimes listed as a/2, this really depends 
on how it is defined. A GaAs unit cell consists of two Ga layers and two A1 layers. 
During growth, it is possible to monitor the Ga and As layers separately. If each of these 
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layers is defined as one monolayer, then there are four monolayers and 1 ML = a/4 by 
definition. On the other hand, if one defines one Ga layer plus one As layer as one 
monolayer, then there are two monolayers and 1 ML = a/2 by definition. The gain is also 
calculated, and the device yields high gain. 




Figure VI.B.2. Results of the Monte Carlo simulation of the 4.4 THz structure, with and 
without including interface roughness scattering, (a) Current density as a function of 
applied electric field bias, (b) Gain as a function of bias. 

This analysis indicates that interface roughness scattering effects are only 
appreciable for roughness greater than ~ 1 ML, and that overestimation in the current 
density may occur assuming a roughness of a/2. The results presented earlier for the step 
well structure, which included interface roughness scattering assuming a roughness of 
a/2, are likely to be overestimating the current density of the structure. While interface 
roughness can vary from sample to sample, this analysis shows the importance of 
maintaining correct tolerances during growth. Because our analysis shows reasonable 
agreement with the experimentally measured current densities, it indicates that our Monte 
Carlo analysis of the step well structures should have similar accuracy. 
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VII. CONCLUSIONS AND DISCUSSION 


This research focused on modeling electron transport properties of THz LO- 
phonon QC structures, including the design and analysis of new step well QC structures. 
From the Monte Carlo simulations, it was found that step well structures are capable of 
high gains and oscillator strengths, high injection efficiencies, with comparable 
characteristics to other square well designs, but do have increased scattering from the 
upper state to the lower sates. The peak gain of ~ 87 cm 1 found in these simulations of a 
step well QC structure, is higher than the peak gain of ~ 73 cm 1 reported for a simulated 
LO-phonon square well design, which also showed similar overestimation of the current 
density from their Monte Carlo simulations. 114 The current density near the injector 
anticrossing was higher in these step well QC structure results and some of the electron 
temperatures were slightly higher as well. This overestimation of the current density was 
likely due to not taking into account coherence effects. 

Perhaps one of the most important results from these simulations on the step well 
QC structures, showed high injection efficiencies are possible using a step well injector. 
Though the injection efficiency was much higher than typical for square well LO-phonon 
designs with similar oscillator strengths (~ 1.6 times higher), the current density was still 
high due to scattering from the upper to lower states. It may be possible to utilize a step 
well injector to improve injection efficiency, even if the step is not used to arrange all 
three energy states within the same well. Utilizing a step well injector with resonant 
tunneling, could yield the benefits of both good injection efficiency and improved upper 
to lower state lifetime. Our analysis is the first to suggest using more than one potential 
height, using two or more A1 compositions, in a QCL active region structure. 

The step well structure that was grown and processed, failed experimentally. 
Because of the large contact resistance and the low current through the device at the 
injector anticrossing, it was determined this was due to a growth or contact fabrication 
problem. Therefore, it was not possible to conclude the practical feasibility of such 
structures. Testing is being conducted to detennine the source of the fabrication problem. 
The step well design remains a promising approach that differs from the two successful 
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approaches that have previously been used (miniband and square well LO-phonon). 
Because the middle state (upper phonon or lower lasing state) is a single energy state, 
contrasting to previous LO-phonon based QCL designs that have doublet states, for long 
wavelength lasers this could be important because it eliminates the possibility for unwanted 
THz absorption that could otherwise occur between those doublet states. 

The Monte Carlo simulations performed on a high power conventional square 
well LO-phonon design, showed reasonable agreement with experimentally measured 
results. This indicated that our simulations of the step well structures, was likely of 
similar accuracy. It was also found that interface roughness scattering becomes significant 
for roughness greater than approximately one monolayer, and including interface roughness 
scattering in simulations can lead to an overestimation of the current density in the device, 
depending on the roughness height chosen. This illustrates the importance of maintaining 
tight growth tolerances for QC structures. 

Although a high oscillator strength can lead to high gain, provided the upper state 
lifetime can be maintained, optimum temperature devices may be realized using lower 
oscillator strength structures. It should be pointed out that the best temperature 
performing devices are due to their higher frequency of operation, and none across the 
THz frequency range perform much better than T ~ E ra diativJkB = Hco/k/;. ’ A longer 
wavelength device might operate at a lower temperature, but with respect to E rac u a ,i Ve /kB be the 
same or better than a shorter wavelength device of higher operating temperature. Even in step 
well structures, it is possible to make the radiative transition more diagonal by utilizing an 
additional well prior to the step, and also by not arranging all three energy levels within the 
same well. 

The overestimation between weakly coupled states that occurs in the Monte Carlo 
simulations, can be dealt with either by including dephasing in a density matrix Monte Carlo 
simulation or by using nonequilibrium Green’s functions. The nonequilibrium Green’s 
function technique is more intuitive from a quantum view point, and has several advantages. 
It can be used to include dephasing, many-body effects, and allow for spectral 
estimation. 115,116 This technique is useful in many other areas of condensed matter physics as 
well. We are continuing to research these techniques and structures discussed, and are also 
looking to include the effects of other parameters such as temperature in our analysis. 
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